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Abstract. This paper presents a new numerical method for the com- 
pressible Navier-Stokes equations governing the flow of an ideal isen- 
tropic gas. To approximate the continuity equation, the method utilizes 
a discontinuous Galerkin discretization on piecewise constants and a 
basic upwind flux. For the momentum equation, the method is a new 
combined discontinuous Galerkin and finite element method approxi- 
mating the velocity in the Crouzeix-Raviart finite element space. While 
the diffusion operator is discretized in a standard fashion, the convection 
and time-derivative are discretized using discontinuous Galerkin on the 
element average velocity and a Lax-Friedrich type flux. Our main result 
is convergence of the method to a global weak solution as discretization 
parameters go to zero. The convergence analysis constitutes a numerical 
version of the existence analysis of Lions and Feireisl. 
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1. Introduction 

In this paper, we will construct a convergent numerical method for the 
compressible Navier-Stokes equations: 

Qt + div{gu) = m(0,T)xQ, (1.1) 
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(gu) t + div(gu <g> u) = - Vp(g) + An, in(0,T)xfi, (1.2) 

where SI C i 3 is an open, bounded, domain with Lipschitz boundary dfl, 
and T > is a given finite final time. The unknowns are the fluid density 
g and vector velocity u, while the operator ® denotes the tensor product 
of two vectors ((v (g) v)ij = ViVj). The mechanism driving the flow is the 
pressure p(g) which is assumed to be that of an ideal isentropic gas (constant 
entropy) ; 

p(q) = CLQ 1 - 

To close the system of equations (1.1) - (1-2), standard no-slip boundary 
condition is assumed 

«|an = 0, (1.3) 

together with initial data 

g(0,-) = g €L"' +1 (Q), gu(0, ■) = m € iJ(O). (1.4) 

From the point of view of applications, the system (1.1) - (1.2) is the 
simplest form of the equations governing the flow of an ideal viscous and 
isentropic gas [1, 22]. In the available engineering literature, the reader 
can find a variaty of flows for which the assumption of constant entropy 
(reversibility) is a reasonable approximation. However, while viscosity in 
(1.2) is modeled by the Laplace operator (An), practical applications will 
make use of a more appropriate stress tensor, the simplest being that of a 
Newtonian fluid with constant coefficients 

div § = \i div (Vu + Vu T ) + AV div u. 

Note that our diffusion is a special case of §. Indeed, S reduces to the 
Laplace operator when /i = 1 and A = ^. The analysis in this paper can be 
generalized to all relevant cases of constant non-vanishing u and A. How- 
ever, this comes at the cost of more terms as the chosen finite element space 
(Crouzeix-Raviart) does not satisfy Korn's inequality [19], but with a negli- 
gible gain in terms of ideas and novelty. The more physically relevant case 

7-1 

where fi and A are functions of the density of the form g 2 is not included 
in the available existence theory (cf. [4, 13]) and there does not seem to be 
any obvious way of including it here either. 

In terms of physical applicability of the forthcoming results, a more press- 
ing issue is the equation of state for the pressure. For purely technical rea- 
sons, we will be forced to require that 7 is greater than the spatial dimension; 

7 > 3. 

This is a severe restriction on 7 with no physical applications (to the au- 
thor's knowledge). Kinetic theory predicts a value of 7 depending on the 
specific gas in question: ~ | for monoatomic gas (e.g helium), ~ | for di- 
atomic gas (e.g air), and creeping towards one for more complex molecules 
and/or higher temperatures. It should be noted that global existence is only 
known for 7 > | and hence only in the case of monoatomic gas (see [10] 
and the references therein). In this paper, the restriction on 7 seems ab- 
solutely necessary to prove convergence of the method, but is not required 
for stability. In fact, the strict condition on 7 is related to the numerical 
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diffusion introduced by the method and is in perfect analogy to the prob- 
lems encountered in a vanishing diffusion limit for the continuity equation 
(1.1) (see for instance [ ]). That being said, it might be calming that the 
upcoming analysis can be easily extended to pressures of the form 



where k > may be chosen arbitrarily small. Hence, for all practical pur- 
poses, the numerical method presented herein is convergent for cases very 
close to the physically relevant ones. 

The numerical literature contains a vast body of methods for compressible 
fluid flows. Many of them are widely applied and constitute an indispens- 
able tool in a variety of disciplines such as engineering, meteorology, or 
astrophysics. While the complexity and range of applicability of numerical 
methods for compressible viscous fluids is increasing, the convergence prop- 
erties of any of these methods have thus far remained unknown. In fact, 
prior to this paper, there have not been reported any convergence results 
for numerical approximations of compressible viscous flow in more than one 
spatial dimension. In one dimension, the available results are due to David 
Hoff and collaborators [5, 26, 27, 28] (see also [13]) and concerns the equa- 
tions posed in Lagrangian variables and relies on the ID existence theory 
which have yet to be extended to multiple dimensions. That being said, in 
the recent years there have appeared a number of convergent methods for 
simplified versions of (1.1)-(1.2). In [8, 9, 11], convergent finite volume and 
finite element methods for the stationary compressible Stokes equations was 
developed. Simultaneously, in [14, 15, 16], K. Karlsen and the author devel- 
oped convergent finite element methods for the non-stationary compressible 
Stokes equations. In the upcoming analysis, we will utilize ideas from all of 
these recent papers. 

Let us now discuss the choice of numerical method. For the approximation 
of the continuity equation, we will use the standard upwind finite volume 
method. However, we will formulate this method as a discontinuous Galerkin 
method where the density g is approximated by piecewise constants. For the 
velocity we will use the Crouzeix-Raviart finite element space. The method 
and some its properties was originally inspired by [25, 18] (cf. [14]), but vari- 
ants can be found several places in the literature (see for instance [12]). For 
the momentum equation, we are leaning on the knowledge gained through 
[8, 9, 11, 14, 15, 16], from which it is evident that proving convergence for 
any numerical method is a hard task. In particular, to develop a numeri- 
cal analogy of the continuous existence theory, it seems necessary that the 
numerical discretization of the diffusion and pressure respects orthogonal 
Hodge decompositions (see Section 8 for an explanation) in some form. 

The distinctively new and completely original feature of the upcoming 
method is the discretization of the material derivative in the momentum 
equation. Our discretization will be of the form 

(gu)t + div(gu (g) u) 



p(g) = ag 11 + Kg 



.72 



7i e (1,2), 7 2 > 3, 
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where v~h denotes the I? projection of Vh onto piecewise constants. The 
operator \Jp(gu <S> u) is the specific upwind flux which we shall use (see 
Section 3 for precision). Hence, the material derivative is discretized using 
discontinuous Galerkin with the same polynomial order as the continuity 
discretization. This stands in contrast to the diffusion and pressure terms 
which are solved with the Crouzeix-Raviart element space and hence in par- 
ticular with piecewise constant divergence div Vh matching the density (and 
pressure) space. In the language of finite differences, this means that the 
pressure and diffusion is solved using staggered grid while the material de- 
rivative are solved at the same nodal values as the density. The great benefit 
of this approach is that it solves the long-lasting problem of incorporating 
both the hyperbolic nature of the material derivative and the nature of the 
diffusion-pressure coupling. In particular, our main result yields stability 
and convergence for all Mach and Reynolds numbers. 

By proving convergence of a numerical method we will in effect also give 
an alternative existence proof for global weak solutions to the equations 
(1 . 1)-(1 .2) . While the first global existence result for incompressible Navier- 
Stokes was achieved by Leray more than 80 years ago, a similar result for 
compressible flow was obtained by P-L. Lions in the mid 90s. In the cele- 
brated book [ ], Lions obtains weak solutions of (1.1)-(1.2) as the a.e weak 
limit of a sequence of approximate solutions. Consequently, from the point 
of view of analysis, we will in this paper perform similar analysis to that 
of [17], but for a numerical method. That is, we will develop a numerical 
analogy of the continuous existence theory. The key difficulty when passing 
to a limit is presented by the non-linear pressure p(g). Specifically, com- 
pactness of g is needed, while the available estimates provides no form of 
continuity of g. In all current proofs of existence, the necessary compactness 
is derived using renormalized solutions together with a remarkable sequen- 
tial continuity result for the quantity p(g) — divu. The result provides a.e 
convergence of the density, but gives no insight into continuity properties 
of g. In the original proof, Lions needed that 7 > §■ The existence theory 
was further developed by Feireisl and the requirement lowered to 7 > |. 
This seems to be optimal in the absence of pointwise bounds on the den- 
sity. However, it is interesting that this still does not include the case of 
air in three dimensions. The reader is strongly encouraged to consult [23] 
for a thorough and well- written introduction to the mathematical theory of 
solutions to (1.1)-(1.2). 

Organization of the paper: In the next section, we will go through some 
preliminary knowledge where we attempt to make clear the solution concept, 
basic compactness results, and the finite element theory used in the analy- 
sis. Then, in Section 3, we present the numerical method, give some basic 
properties of the method, and state the main convergence result (Theorem 
3.5). We will then move on to Section 4 in which we establish stability of the 
method and draw conclusions in terms of uniform integrability. The follow- 
ing section, Section 5, is a fundamental section where we will estimate the 
weak error (in a weak norm) of the transport operators in the discretiza- 
tion. The material contained in this section will be used ubiquitously in 
the convergence analysis. In Section 6 we prove higher integrability of the 
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density. That is, the density enjoys more integrability than the energy esti- 
mate provides. Then, in Section 7, we will pass to the limit in the method 
and conclude that the limit is almost a global weak solution. It will then 
only remain to prove strong convergence of the density approximation. In 
Section 9, we establish the fundamental ingredient in the proof of density 
compactness; - weak sequential continuity of the effective viscous flux. Fi- 
nally, in Section 10 we prove strong convergence of the density and conclude 
the main result (Theorem 3.5). The paper ends with an appendix section 
containing the proof of well-posedness for the numerical method. 

2. Preliminary material 

The purpose of this section is to state some results that will be needed in 
the upcoming convergence analysis. 

2.1. Weak and renormalized solutions. 

Definition 2.1. We say that a pair (g,u) is a weak solution of (1.1) - (1-2) 
with initial condition (1.4) and boundary condition (1.3) provided: 

(1) The continuity equation (1.1) holds in the sense of distributions 

/ f?(<£t + u ' V0) dxdt = — I Qo4>(0, •) dx, 
/o Jn Jn 

for all (f> G C§°([0,T) x fi). 

(2) The momentum equation (1.2) holds in the sense 

/ / —Quv t — gu® u : Vv — p(g) div v + VuVv dxdt 
Jo Jn 

m,QV (0, •) dx, 

n 

for all v G [Cg°([0,T) x ft)] 3 . 

(3) The energy inequality holds 

,2 



f Q U 1 t \ r 

sup / — — I -p{Q) dx 

e(o,t) Jn 1 7 - 1 

+ / / |Vn| 2 dxdt < I ^ + -^—p(go) dx. 
Jo Jn Jn 2 Qo 7-1 



Definition 2.2 (Renormalized solutions). Given u G L 2 (0, T; W Q ' 2 (ft)), we 
say that g G L°°(0, T; L 7 (ft)) is a renormalized solution of (1.1) if 

B(g) t + div (B(g)u) + b(g) div u = 0, 

inV ([0,T) x 0) for any B G C[0, oo)nC 1 (0, oo) with 5(0) = and b(g) := 
B'(g)g-B(g). 

We shall need the following well-known lemma [17] stating that square- 
integrable weak solutions g are also renormalized solutions. 

Lemma 2.3. Suppose (g, u) is a weak solution according to Definition 2.1. If 
g G L 2 ((0, T) x Q)), then g is a renormalized solution according to Definition 
2.2. 
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2.2. Compactness results. In the analysis, we shall need a number of 
compactness results. 

Lemma 2.4. Let X be a separable Banach space, and suppose v n : [0, T] — > 
X*, n = 1,2, ... , is a sequence for which ||^n||_£,o°([o t]-x*) — ^> for some con- 
stant C independent ofn. Suppose the sequence [0,T] 3 t h- > {v n (t), $)x*,X> 
n = 1,2, . . . , is equi- continuous for every <1> that belongs to a dense subset 
of X. Then v n belongs to C ([0, T]; X* eak ) for every n, and there exists a 
function v G C ([0, T]; -^ eak ) swc/i that along a subsequence as n — ^ oo t/iere 
holds v n ^v inC([0,T];X* eak ). 

To obtain strong compactness of the density approximation, we will utilize 
the following lemma. 

Lemma 2.5. Let O be a bounded open subset of ~R M , M > 1. Suppose 
g: M — t- (—00,00] is a lower semicontinuous convex function and {v n } n>1 is 
a sequence of functions on O for which v n ^ v in L 1 (0), <7(f n ) E -^ 1 (0) /or 
eoc/i n, g(v n ) — 1 5(f) in L 1 (0). T/ien 5(f) < g(i>) a.e. on O, g{v) G L 1 (0), 
and J G <?(?;) dy < liminf n _ 5>00 j g{v n ) dy. If, in addition, g is strictly convex 

on an open interval (a, b) C M and g(v) = g(v) a.e. on O, then, passing to a 
subsequence if necessary, v n (y) — > v(y) for a.e. y £ {y £ O \ v{y) £ (a, 6)}. 

In what follows, we will often obtain a priori estimates for a sequence 
{v n } n>1 that we write as "v n Gb X" for some functional space X. What 
this really means is that we have a bound on H^nllx that is independent of 
n. The following lemma is taken from [14]. 

Lemma 2.6. Given T > and a small number h > 0, write (0, T] = 
Uj^fa-i.t*] with t k = hk and Mh = T. Let {f h }^ >0 , {g h }h>o be two 
sequences such that: 

(1) the mappings t 1— > g^{t,x) and t \- > fh(t,x) are constant on each 
interval (tk-ijtk], k = l,...,M. 

(2) { fh}h>o and {g h }h>o converge weakly to f and g in L Pl (0, T; L qi (O)) 
and L P2 (0, T; L q2 (O,)), respectively, where 1 <pi,qi < 00 and 

1111 

— + — = - + - = 1. 

Pi P2 qi q2 

(3) the discrete time derivative satisfies 

9h(t,x)-g h (t-h,x) £b Ll w -l,i m 
h 

(4) \\fh{t,x)- f h (t,x-^)\\ LP2{ Q yT . Lq2lsl)) as |£| 0, uniformly in h. 
Then, gnfh gf in the sense of distributions on (0,T) x £1. 

2.3. Finite element spaces and some basic properties. Let Eh denote 
a shape regular tetrahedral mesh of Vt. Let Th denote the set of faces in Eh- 
We will approximate the density in the space of piecewise constants on Eh 
and denote this space by Qh(Q). For the approximation of the velocity we 
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use the Crouzeix-Raviart element space [6]: 

V h (n) = <L G L 2 (Q) : v h \ E G P?(E), V£ G £7 h , 



(2.2) 



Kir <*s(x) = o, vr g r fc 



where [-| r denotes the jump across a face T. To incorporate the boundary 
condition, we let the degrees of freedom of V^(O) vanish at the boundary. 
Consequently, the finite element method is nonconforming in the sense that 
the velocity approximation space is not a subspace of the corresponding 
continuous space, W ' 2 (fi). 

For the purpose of analysis, we shall also need the div-conforming Nedelec 
finite element space of first order and kind [20, 21] 



tf h (Sl) = { v h , dfrv h G L 2 (n) : v h \ E G Pq ® F^x, VE G E h , 

[v ■ uj dS(x) = 0, Vr G T h 



L 



(2.3) 



'r 

We introduce the canonical interpolation operators 

u v h : wl' 2 (p) ^ v h (n), 



n« : L 2 (fl) i ^ Q h (Sl), (2.4) 



defined by 



J Rh v h dS(x) = J v h ds(x), vr g r h) 

J ll%v h -udS(x) = J v h -udS(x), VTeT h , (2.5) 

/ n®(f) dx= I <j) dx, \/E G E h . 
Je Je 

Then, by virtue of (2.5) and Stokes' theorem, 

divD^v = divfclljfw = Iljjdivi;, cml h Ii\v = njj curlw, (2.6) 

for all v G VFg' 2 (r2). Here, curl/j and div^ denote the curl and divergence 
operators, respectively, taken inside each element. 

Let us now state some basic properties of the finite element spaces. We 
start by recalling from [3, 6, 20] a few interpolation error estimates. 

Lemma 2.7. There exists a constant C > 0, depending only on the shape 
regularity of Eh, such that for any 1 < p < oo, 

Lv{E)i 

linjfu - v\\ LP(E) + h\\V(Ulv - v)\\ LP(E) < ch s \\V s v\\ LP{E) , s = 1, 2, 
\\Iltfv-v\\ LP{E) +h\\div(ltffv-v)\\ LP{E) < ch s \\V s v\\ LP{E) , s = 1,2, 
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for all <j> e W 1 ' P (E) and v £ W S ' P (E). 

By scaling arguments, the trace theorem, and the Poincare inequality, we 
obtain 

Lemma 2.8. For any E € Eh and <ft £ W 1,2 (E), we have 

(1) trace inequality, 

UWlhd < ch~* (UWlhe) + h E \\V<l>\\v{E)) , vrer.n dE. 

(2) Poincare inequality, 
-tL I <t>dx <Ch E \\V^\\ L2(E) . 

1-^1 JE L 2 (E) 

In both estimates, He is the diameter of the element E. 

Scaling arguments and the equivalence of finite dimensional norms yields 
the classical inverse estimate (cf. [2]): 

Lemma 2.9. There exists a positive constant C, depending only on the 
shape regularity of Eh, such that for 1 < q,p < oo and r = 0, 1, 

\\<Ph\\ W r, P{E) < c/r r+mm{0 'p _ 5 } Uhh^E) > 

for any E S Eh and all polynomial functions <f>h G Fk(E), k = 0, 1, . . .. 

Since the Crouzeix-Raviart element space is not a subspace of W 1,2 , it 
is not a priori clear that functions in this space are compact in LP , p < 6. 
However, from the Sobolev inequality on each element and an interpolation 
argument we obtain the needed result. 

Lemma 2.10. For Uh G Vh(£l), let p satisfy 2 < p < 6 and determine a 
such that 

1 a (l-o) 

- = - + -. 

p 2 6 

The following space translation estimate holds 

a 

\\u h (-) ~ U h {- - OllLP(n\{a«dist(a!,«l)}<[f|) < C { h + ICl ) 2 || Vtt fo || Z 2 (a ) , 

where C is independent of h and £. 

Proof. From the work of Stummel [24] , we have that 

IMO - u h {- - 011^(0) < C + l|V^|| L2( n). (2.7) 
The standard Sobolev inequality gives 

IMO - u h {- - Oh^n) < C\\V h u h \\mn)- (2-8) 
Hence, the proof follows from interpolation between (2.7) and (2.8). 

□ 

Finally, we recall the following well-known property of the Crouzeix- 
Raviart element space. 

1 2 

Lemma 2.11. For any Uh 6 Vft(fi) and v £ W ' (fl), 

[ V h u h V h (ulv-v) dx = 0. (2.9) 
Jn 
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Proof. By direct calculation, we obtain 



Vh^V /l (n^v - v) dx 



n 




[ Au h (U^v - v) dx + [ (Vu h ■ u){n\v - v) dS(x) 



£ JE JdE 



Now, since Uh is linear on each element Auh = and Vu^ is constant. 
Moreover, since the normal vector v is constant on each face of the element, 
we have that 



by definition of the interpolation operator IlY. Hence, we have proved (2.9). 



For a given timestep At > 0, we divide the time interval [0, T] in terms 
of the points t k = kAt, ft = 0, ... , M, where we assume MAt = T. To dis- 
cretize space, we let {-E^l/i be a shape-regular family of tetrahedral meshes 
of O, where h is the maximal diameter. It will be a standing assumption 
that h and At are related like At = ch for some constant c. We also let 
denote the set of faces in E^. 

The functions that are piecewise constant with respect to the elements of 
a mesh are denoted by Qh(ty and by V/i(f2) we denote the Crouzeix- 
Raviart finite element space (2.2) formed on E^. To incorporate the bound- 
ary condition, we let the degrees of freedom of V/j(0) vanish at the boundary: 



We will need some additional notation to accommodate discontinuous 
Galerkin discretization. Related to the boundary dE of an element E, we 
write /+ for the trace of the function / taken within the element E, and /_ 
for the trace of / from the outside. Related to a face T shared between two 
elements E- and E+, we will write /+ for the trace of / within E + , and 
/_ for the trace of / within EL. Here E- and I?+ are defined such that v 
points from E- to E+, where v is fixed as one of the two possible normal 
components on T. The jump of / across the face T is denoted lf} r = /+—/_. 

To pose the method, and in the convergence analysis, we shall need the 
canonical interpolation operators (2.4). In fact, we shall need the operators 
and to such an extent that we introduce the convenient notation 




□ 



3. Numerical method and main result 




v = ni% 



(3.1) 



To discretize the convective operator div(^u) in the continuity equation 
(1.1), we will utilize a standard upwind method in the degrees of freedom 
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Eh - the mesh. 

E - an element in the mesh. 

dE - the boundary of E. 

r - a face in the mesh. 

Th - all faces in the mesh. 

Qh(Q) - the space of piece constant scalars on E^. 

Vh(fl) - the Crouzeix-Raviart vector space on E^. 

Nh{ty - the div conforming Nedelec space 

of first order and kind. 

11^ - the I? projection operator onto Qh- 

11^ - the canonical interpolation operator onto Vh- 

11^ - the canonical interpolation operator onto Mh- 

f - n^/ (the piecewise constant projection). 

v - IT^ v (the Nedelec interpolation) . 

/+ - max{/,0}. 

/- - min{/,0}. 

f+\dE - the trace of / taken from within E. 

f-\dE - the trace of / taken from outside E. 

/+|r - the trace of / taken against the normal vector v. 

/— |r - the trace of / taken in the direction of v. 

Mr - /+-/- 

Up(Qu)\eE - Q+(u h - V) + + Q-(u h - V) . 

Up(^u)| r - g-(u h ■ v) + + g+(u h ■ v)~ . 

XJp(gu<giu)dE - Up + (^n)n + + XJp~(gu)u-. 

XJp(gu<giu)r - Up - (gu)u + + XJp + (gu)u-. 

Table 1. Notation 



of Uh- The upwind discretization is defined as follows 

Up(£»n)| r = Q- ^j^i j^u h ■ v dS(x)j + g+ ^-^ J^u h ■ v dS(x) j , 

= g-(u h ■ u) + + Q+(u h ■ v)~ , VTeT h , 

where we have used the notation (3.1). 

For the convective operator div(gu (g) u) in the momentum equation (1.2) 
we will use the following Lax-Friedrich type upwind flux 

\Jp(gu ® u) = Up + (gu)u+ + Up~ (gu)u- . (3-3) 

Observe that this operator is posed for the average value of over each 
element. This is non-standard and has to the author's knowledge not been 
studied previously. We are now ready to pose the method. 

Definition 3.1 (Numerical method). Let go E L 7+1 (0) and mo E La (J)) be 
given initial data and assume that T > is a given finite final time. Define 
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the numerical initial data 



el 



7/ 

u h 



n 



V 



m 



go + nh 

where k is a small positive number. Determine sequentially 



(3.4) 



(g k h ,u k h )eQ h (n)xV h (n), 
satisfying, for all q^ G Qh^l), 



k = l, 



,M, 



I 

JQ 



si 



At 



q h dx-Y^ / Up fc (£u) lq h j r dS(x) 



+ /ll ~ £ E/[^] r Mr dS(x) = 0, 



(3.5) 



and for all Vh £ Vh(Q), 



At 



r 



+ 



Vh dx 

r 

V h u k h V h v h - p(g k h ) &\vv h dx 

i 



Up k {gu®u) lv} r dS(x) 



(3.6) 



lv h j T dS(x) = 0, 



where e > g should be chosen as small as possible. 

For the purpose of analysis, we will need to extend the pointwise-in-time 
numerical solution (^,u^), k = 0, . . . ,M, to a piecewise constant in time. 
For this purpose, we will use the following definition 

(Q h ,u h )(t,-) = (Q k h ,u k h ), te[t k ,t k+1 ), k = Q,...,M. (3.7) 

Remark 3.2. The terms involving h l ~ e in (3.5) - (3.6) are needed for purely 
technical reasons. In particular, they are not needed to obtain stability of 
the method or to prove convergence of the continuity method (3.5). See 
Section 8 and [9] for more on why they are needed. 

3.1. The method is well-defined. Since the numerical method is nonlin- 
ear and implicit, it is not trivial that it is actually well-defined. In addition, 
the transport operators in the momentum equation is posed for the element 
average velocity and hence does not provide a full set of equations in them- 
selves. In fact, it is only due to the diffusion that the system has a sufficient 
number of equations for the degrees of freedom of Uh- We shall prove the 
existence of a numerical solution through a topological degree argument. 
The proof is very similar to that of [12, 14] and is for this reason deferred 
to the appendix. 



Proposition 3.3. For each fixed h > 0, there exists a solution 

(g k h ,u k h )eQ h (n)xv h (n), e k h (-)>o, k = i,...,M, 

to the numerical method posed in Definition 3.1. 



In the upcoming analysis, we will need that the density solution is positive. 
For this purpose, we recall the following result from [14] (see also [12]): 
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Lemma 3.4. Fix any k = 1, ...,M and suppose g k h 1 G Qh(ty> u \ ^ 
Vh(Cl) are given bounded functions. Then the solution g\ 6 Qh(^) o/ the 
discontinuous Galerkin method (3.5) satisfies 



mm g h > mm o. 



1 



x&n h x&i h \ l + Ai|| divfeu^lioo^) / ' 
Consequently, if g^~ l (-) > ; i/ien £?|(-) > 0. 

3.2. Main result. Our main result is that the numerical method converges 
to a weak solution of the compressible Navier-Stokes equations (1.1) - (1-4). 

Theorem 3.5. Suppose j>3,T>0isa given finite final time, and that 
the initial data (go, mo) satisfies 

/ + Mgo) dx <C, go £ L°°(0,T;L^ +1 (Cl)). 

Jn 2^0 7-1 

Let {(ghiUh)}h>o be a sequence of numerical solutions constructed according 
to Definition 3.1 and (3.7) with At = ch. Along a subsequence as h — > 0, 

u h ^u in L 2 (0,T;L 6 (Cl)), 

V h u h Vu in L 2 (0, T; L 2 (Cl)) 

QhUh, QhUh QU in L (0, T; L^+i (Cl)), 

QhUh QU in C(0, T; L^+i (CI)). 

QhUh ®Uh^gu®uinL (0, T; L 3 +t (Cl)), 

Qh g a.e and in L p loc ((0,T) x Cl), p < 7 + 1, 

where (g, u) is a weak solution of the isentropic compressible Navier-Stokes 
equations (1.1) - (1-2) in the sense of Definition 2.1. 

The proof of Theorem 3.5 will be developed in the remaining sections of 
this paper. The final conclusion will come in Section 10.1. 

4. Energy and stability 

In this section we will prove that our method is stable and satisfies a nu- 
merical analogy of the continuous energy inequality (2.1). Both the stability 
estimate and large parts of the subsequent convergence analysis relies on a 
renormalized type identity derived from the continuity scheme (3.5). The 
proof of this identity can be found in [14, 15]. We shall utilize the following 
form: 

Lemma 4.1. Let (g h ,u h ) solve the continuity scheme (3.5). Then, for any 
B £ C 2 (R+) with B" > 0, there holds 



I 

Jn 



B(Q k h ) dx + (At) -i J B „ {e%){Q k _ gk -i ? dx 



n 

2 



+ 2J/ S/, G*)HJ (KM + fc 1 - 6 ) ds(x) 

< - / (qB'(q) - B(g))divu h dx, 
Jn 
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where g± and g* are some numbers in the range [g k ~ 1 ,Q k ] and [Q-,Q+], 
respectively. 

We now prove our main stability result. 

Proposition 4.2. For given At, h > 0, let (u^, gfy, k = 0, ...,M be the 

numerical approximation of (1.1)-(1.2) in the sense of Definition 3.1. Then, 

n m I C~m 1 2 i 

6hl hl +^-r P ( T)dx 



max 

m=l,...,M Jn 



7-1 



M . 5 M 

+ AtJ2 / \V h u k h \ 2 dx + AtJ2Y, D i ( 41 

k=l Jn i=l fe=l 



< 



I —0 1 2 
8h\ u h\ 



1 



.-> + 7P(e°) dx > 

z 7 — 1 

where the numerical diffusion terms are given by 



n 



^ = E / RpV) 



\u\ ■ v\ dS(x), 



dS(x), 



D\ 
Dl 



^"E r l PttWr dS(x)dt, 
p Jo Jr 

(At)" 1 f P"{Q X ){e k h -Q k H l ? dx, 
Jn 

(Aty 1 [ 
Jn 



el V*-"* 'I 2 dx. 



Proof. Let Vh = u\ in the momentum scheme (3.6), to obtain 



k-l^k-l 



At 



ul dx+ I \V h u k \ 2 dx 



y2 / Up fc (£»n (8) u)u k + dS{x) + / p{g k h ) div u k h dx 
v JdE Jn 



(4.2) 



From Lemma 4.1 with -B(z) = we see that the pressure term 

p(Qh) divu h dx 



> 



(4.3) 



dx-D k -D%- D k . 



7-1 At 

Next, we turn our attention to the first term after the equality in (4.2). 
From the definition of Vp(gu (g) u), we have that 



y2 Up k (gu®u)u k l dS(x) 

w JdE 
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V / (Up+ \gu)u k + + \Jp-(gu)u k _) u k + dS{x) 

^ JdE V ' ' 

JdE l 



<k \2 



2 

fc -fc 



+ Up + (gn) ^ + Up (£>u)2_ 2^5(2;). 

By setting ^ = (l/2)(u k i ) 2 in the continuity scheme (3.5), we see that the 
first term after the equality in (4.4) appears 



k _ fc-l \gk\2 

2/i ^/i \ u h\ 



dx 



At 2 



k \2 



8E 



E / u p fc (^) > 



2 

k \2 







1 


Qh] 


JdE 





k \2 



dS(x) 



dS{x) 



(4.5) 



^+±^ ] u\ dS(x). 



To see the contribution of the second term in (4.4), we first recall that 

Up + (gu)\ dE+ = - Vp-(gu)\ dE _ , dE + n 8E~ = T, 

since the normal vector has opposite signs. Using this, we obtain 

(v k ) 2 (v k ) 2 

- Up-(Qu)^- + Up- (gu)u k _u k + dS(x) 



E / u p + (^) 



(u k ) 2 _. .(u^) 2 tt Ju k + ) 



k \2 



Up (qu)- 



(u k ) 2 

Up + (gu)^^- + \Jp-(gu)u k _u k + - XJp + (gu)u k + u k _ dS{x) 



E/iVmi 



(ui) 2 + (u^) 2 



dS{x) = Dl 



2 2 

By applying this together with (4.5) in (4.4), we discover 
~E/ J Jp k (gu®u)u k r dS(x) 

E JdE 

+ ^E 



Qh 

8E 

fc-l 10*12 



) ^ + dS(x) 



(4.6) 



At 2 2 
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Now, setting (4.6) and (4.3) into (4.2) reveals 

ft%-tfr l tfr\ k _ ej - eV \uj? i P(fi) - Piet 1 ) dx 
n At h At 2 7-1 At 

+ f \V h u k h \ 2 dx + D\ + D\ + D\ + D\ < 0. 
Jn 

A simple calculation gives 



/ 



kc:k _ h-X-rjk-l k _ fc-1 \Cjk\2 

°h u h 6h U h u k _ °h Sh \ u h\ j 

At h At 2 

Jci fe|2 _ Jfe-l| fc-l|2 _ 7 , fc - 1 | 2 

2At At 
fci fe|2 _ k-l\ u h-Ul 

^ A h K 1 dx + Z>*. 

2At 5 

Consequently, by combining the two previous inequalities, 

eiWjf - e k h - l \u k h - l \ 2 M) - KgT) 

2At (7-1) At 



5 

/ |v,4l 2 dx + ^A fc <o. 



We conclude by multiplying with At and summing over k = 1, . . . , M. □ 

Observe that the energy estimate does not provide L°° control in time 
on Qh\uh\ 2 - Instead, we only gain this control on the projection Qh\uh\ 2 - 
The following corollary is an immediate consequence of the energy estimate 
(Proposition 4.2) and the Holder inequality (cf. [10]). 

Corollary 4.3. Under the conditions of Proposition 4-2, 

8h L°°(0, T; L 7 (il)), p(g h ) e b L°°(0, T; L 1 ^)) 
u h € b L 2 (0, T; L 6 (ft)), V h u h G 6 L 2 (0, T; L 2 (fl)), 

QhU h e b L°°(0,T;L m °°(f})), g h u h e h L 2 (0, T; L m2 (VL)), 
Qh \u h \ 2 SL°°{Q,T;L\VL)), Q h \u h \ 2 G 6 L 2 (0,T;L C2 (ft)), 
where the exponents are given by ( since 7 > 3) 

27 3 67 37 3 

™°o = — tt > 7T> ™2 = — — > 3, c 2 = — — > -. 
7+12 3+7 3+72 

5. Estimates on the numerical operators 

To prove convergence of the numerical method, our strategy will be to 
adapt the continuous existence theory to the numerical setting. We will 
succeed with this by controlling the weak error of the numerical operators 
relative to their continuous counterparts. The purpose of this section, is to 
derive the needed error estimates. 

For notational convenience, let us define 

/(t)-/(t-At) 
Ut! ~ At 
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and observe that this satisfies 



DtQh(t) 



5.1. The convective discretizations. We begin by deriving identities for 
the distributional error of the numerical convection operators. 

Lemma 5.1. Fix two functions (ft 6 C^°(Q), v G [Cg°(Jl)] <i . The numerical 
transport operators in (3.5) and (3.6) satisfies the following identities 



£jf Up(gu) 

/ Up(^« ® s) 





d5(x) = f 


r 




d5(x) = f 


r ./n 



(5.1) 



8=2 



(5.2) 



where the error functionals Pi, i = 1, ... ,4, are 

Pi{v) Qh divu h u h (H%v - v) dx. 

E Je 

Proof. Using the continuity of Vp(gu) and (ft across edges, we calculate 



E/ u pM 



n 



dS(x) 



V / Vp(gu) dS(x) 

^ JdE 

V / Vp(gu) Utfcft-A dS(x) 
E JdE v ' 

I ■ v) + + g-(u h ■ i/)-) (njfy - </>) ^(x), 
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where the last identity is the definition of \Jp(gu). Next, we add and subtract 
to deduce 

^jTupCH [njfy] r dS(x) 

+ ( e _ - Q+ )(u h ■ v)- (njfy - </>) dS(x) 

= y QhUhS?<t> dx - ^2 J Sh divu h (lL®<f> - (pj dx + Pi(<j>). 



We conclude (5.1) by recalling that div/j Uh is constant on each element and 
hence the second term is zero. 

To derive (5.2), we apply the definition (3.3) of \Jp(gu ® u) and add to 
subtract to obtain 



QU ® U) 



dS{x) 



J \Jp(gu®u) (nlv- v) dS(x) 

E 

V / (Up + (eu)u + + Up_(e«)«_) fn^-u) d5(x) 

^ I (XJj>{qu)u + + Up_(^n)(n_ -«+)) (nj^ - v) dS(x) 

^1 ■/ dE 

V / Up(^n)S + (LT> - v) dS(x) + P 3 (v). 

^7 JdE v 7 



To proceed, we apply the definition of \Jp(gu) (3.2) and add and subtract 
to obtain 



/ Up(£u ® 



dS(x) 



V / {Q+(u h ■ v) + + g^(u h ■ 1/)-) u+ (itfv - u) d5(s) + P 3 (* 

V / g+(u h ■ v)u+ [H\v - v) 

+ (0- ~ Q+)(u h ■ v)~u+ (uXv - vj dS(x) + P 3 (v) 
[ Q+(u h • v)u + (nf ; - v) dS(x) + P 2 (v) + P 3 (v) 
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We can then apply the divergence theorem to the first term to obtain 



<X> u) 



dS(x) 

r 



^2 J^Vp(gu< 

= / QhUh ® uh : Vv dx - / Sh div u h u h (ji\v - v) dx 
Jn ^ J e 

+ P 2 (v) + P 3 (v) 

= / Q h u h ®u h :X7v dx + P 4 (v) + P 2 {v)+P 3 (v), 
Jn 

which is (5.2). □ 



Identities like (5.1) and (5.2) can be derived for any numerical method. 
The difficult part is to control the error terms, in our case is given by Pi, i = 
1, ... ,4. The following proposition provides sufficient control on the error 
terms. Note in particular the integrability required of the test-functions. 

Proposition 5.2. Let (^uf), k = 1,...,M be the numerical solution 
obtained using the scheme (3.5) - (3.6). Let (gh,Uh) be the piecewise constant 
extension of (g^ufy, k = 1,...,M in time to all of [0,MAi] (i.e (3.7)j. 
Then, the Pi, % = 1, ... ,4 in Proposition 5.1 are also piecewise constant in 
time and there exists a constant C > 0, independent of h and At, such that 



Vi(0)l^<^- mm{3 ^'° } C||V^|| L4(OiT;i¥(Q) , (5.3) 

T 

\P 2 {v)\ dt < h^ {3 ^r' 0} C\\Vv\\ L oo i0iT , L3m , (5.4) 

T 

\P 3 (v)\ dt < h*C\\Vv\\ 6 7 , (5.5) 

L4(0,T;L5^(Q)) 
T 2(7-9/4) 

\Pa{v)\ dt<h -i C[|Vu|| iO o( 0) T;£,7(n))- (5.6) 





Ln particular, for 7 > 3, we have that 







T 







\Pi{6)\ dt < h*C\\V6\\ „ 

1 — 11 ^ ll L 4 (0,T;LS-(ny 



/ \P 2 {v)\ + \P 3 {v)\ dt < Chi\\Vv\\ L oo { Q )T . Lim . 
Jo 



(5.7) 



Proof. We will prove one inequality at the time. 
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1. Bound on P£: An application of the Cauchy-Schwartz inequality yields 



f 

Jo 



\Pi\ dt:= I W { 6h \ dE (u h ■ u)-(U^<j> - cf>) dS( 3 

JO £ JdE 

<([ ^2 [ {QhfdE K ■ v \ dS(x)dt ] 

\Jo £ JdE J 

x ( r e i \ u h\ - ^ 

\Jo p JdE 



dt 



(5.8) 



h x 

By setting B{z) = in Lemma 4.1 and integrating in time, we obtain 

h< \ \ (Qh) 2 dx+ I I (g h ) 2 div h u h dxdt 
1 Jn Jo Jn 

< C||£ , /iIIl=o(0,T;L 4 (C)) + ^H^IlL°°(0,T;L 4 (n)ll ^ v h Uh\\ L 2 (0,T;L 2 (Q)) 

< C|l£'/i|lL°°(0,T;L 4 (n))( 1 + ^T)- 

By applying the Trace Lemma 2.8 and the Holder inequality, we deduce 

fT 



h<h 



-i 



< br x C 



o Jn 

T 



|u ft |(nj^- <t>Y dxdt 

\ u h\\LS(U) 

x ( \\H$6- 



(5.10) 



\\ 2 +/l 2 ||V0||\ la \dt 



<^C||^IU 2 (0,T ; La(n))l|V^||^ (0)T;i¥(n))) 

where the last inequality is the error estimate on the interpolation error of 
LT^ in LP (see Lemma 2.7). 

Consequently, by applying (5.9) and (5.10) in (5.8), we see that 



(5.11) 



\P X \ dt < h*C\\Q h \\ L0O (pF. L 4(a i )\\V<l>\\ L 4(0)T;L ^ (n)) - 
From the standard inverse inequality (Lemma 2.9), we have that 



l^/i||L°°(0,T;L 4 (a)) 



< ^ma x {-3^,0} c 



l^felk°°(o,r;iT(n))) 



(5.12) 



where the last inequality follows from 7 > 3. Applying this in (5.11) finally 
yields 

/Vl dt < h^^c\\ Qh \\^ 

Since Corollary 4.3 provides Qh Gb £°°(0, T; L 7 (Q)), we can conclude (5.3). 
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2. Bound on P 2 : An application of the Cauchy-Schwartz inequality yields 



\P 2 \ dt := / S~) / lQhJ aE (uh ■ v) u h (Il%v -v) dS(x) dt 
JO e JdE 

<[[ JZ [ lQhf 9E \ u h-v\ dS{x)dt ) 

\Jo ^Joe J (5.13) 

x (£^J^u h \\u h \QlXv-v) 2 dS(x)dt) 



■ Vh x v/J 2 . 



From (5.9) and (5.12), we have that 



r imax{ — 3— r^.Ol/^ii 



lL°°(0,T;LT(n))- 



(5.14) 



Using the Trace inequality in Lemma 2.8, the Holder inequality, the stability 
of LT^ in L 6 , and the interpolation error estimate for II in L 3 (Lemma 2.7), 
we deduce 



T 



h<h-/ y m 



dx 



x (l|n^-w||i 3(n) + ft 2 C||Vi;||i3 (n) ) '// 



< br x c 



\ u h\\L^(Sl)\\ u h\\Lfi(Sl) 



(||n^-u||i3 (n) + ^ 2 q|v w ||i 3(f2) 



< /iC||u/ l |||a(o,T;X6(n))ll^' ; ll£«'(0,T;L3(n))) 



where the term involving is bounded by Corollary 4.3. 
Now, setting (5.14) and (5.15) in (5.13) yields 



(5.15) 



\P 2 \ dt < h^- min{3 ^' 0} C\\Vv\\ L ^ T . LHn)) , 



which concludes our proof of (5.4). 
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3. Bound on P% : An application of the Cauchy-Schwartz inequality yields 

-T 

\P 3 (v)\ dt 

E/ u P-MI^la^( n ? n ^-«) dS(x) 

rp J dE 











dt 



< 



T 



JT 



|Up(£m)l IM dS{x)dt 



E 



J dE 



|Up(gu)| (IEfn^«-u) dS(s)dt 



:= ii x I 2 . 

From the energy estimate (Proposition 4.2), we have that 

/ M \ 2 



h = E At£> 2 - c ' 



and hence it only remains to bound i 2 . 

Using the definition (3.2) of Vp(gu) and the Holder inequality 



T 2 



E 

E 



JdE 
T 



|Up(gu)| (n^n^u-u) d5(x)dt 



< 



E 



JOS 



(q+ + 8 -)\u h ■ v\ ( [n£n£>] + - u ) ^(x)dt 







L^{dE) L$+i(dE) 



Ii Q h Ii v h v] - V 



e~t dt. 

LTr^idE) 



To proceed, we apply the trace Lemma 2.8 followed by the inverse estimate 
(Lemma 2.9) 

if < 2h- l C\\ eh u h \\ jn_ 



67 



L 4 (o,r ; L57^(n)) 



+ /i 2 ||Vd| 2 



67 



L 4 (0,T;L5t=B(O)) 



< 2/l CH^ftUftll 6 7 

L 2 (0,T;L B^R(n)) 



67 



L 4 (0,T;L5^B(n)) 



IrrF II 2 
lit U — v\\ 67 

1 ft m L 4 (0,T;L5t=B(O)) 



+ /i 2 ||Vu|| 2 



67 



L 4 (0,T;L 5 T=f J (n)) 



< /lC4||^« ft || 67 ||Vw|| 2 67 , 

L 2 (0,T;LB+7(Q)) L 4 (0,T;L5^(n)) 
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where we in the last inequality have used both the interpolation error of LT^ 
and 11^ (see Lemma 2.7). In particular, we have used the following estimate 

2 



< h 2 C\\Vv\\ 2 



67 

Now, from Corollary 4.3 we have that II £>hii/i II 67 < C, and hence 

L 2 (0,T;L^ (O)) 

we can actually conclude that 



/ IP3WI t&< /»3C||Vu| 
J 



L 4 (0,T;L 5 ^ 13 (n))' 



which is (5.5). 



^. Bound on P^: By direct calculation using the Holder inequality 

T 

\P±{v)\ dt 







< || divftitftHia^T;^^)) fj£ j^\QhUh\ 2 (J^v - v) 2 dxdt^j 

<hC\\QhUh\\ ^7 ||Vf ||l/x>(oT;£7(fi)), 

L 2 (0,T;L^(ft)) ' ' V y; 

where we have used that div^«^ Gb £ 2 (0, T; L 2 (fi)) and the interpolation 
error Lemma 2.7. Now, we apply the inverse estimate in Lemma 2.9 to 
obtain 

- 2(7 ~ 9/4) ^ 

^ WQhUhW ^ <h t C\\Q h u h \\ 67 

L 2 (0,T;L^(C)) L 2 (0,T;L3+^(ft)) 

We conclude (5.6) by recalling that the last term is bounded by Corollary 
4.3. □ 

5.2. The material momentum transport operator. In our proof of 
convergence we will need an identity like (5.2) for the discretization of both 
terms in the discrete material transport operator ((gu)t + div(gu ® it)) com- 
bined. To obtain the desired estimates, we will rely on the following weak 
time-continuity result: 

Lemma 5.3. Let (Qh,Uh) satisfy the energy estimate (4.1). Then, 
||^-(^)(-At)|| if{At)T;i|(n)) <(At)^ 



As a consequence, 



|A(^)|| Lf(AtiT;L 3 (c)) <(At)Vc. 



Proof. From the energy estimate (4.1), we have that 

E / P"{ex)\Q k h -e k h - l ]dx<c, 
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where P"{gt) is determined as the remainder in a Taylor expansion and the 
mean value theorem. In particular, a simple calculation gives 



P"(QtM ~ Q k H l ? dx > 1/(7) / [e k h - Q k h - l Y dx, 
where 1/(7) only depend on 7. Hence, 

E \\Sh ~ 6h 1 !ll7(n) - C - 



(5.16) 



By adding and subtracting, we have that 

11 k^k fc-l~fc-l|| 

\\Q h u h -Q h u h || £S 



\Qh y u h~ u h \+ u h\Qh-Qh Jiy (n) 



fc— 1 1 



- Mi *H£3(n) 

C||%lk 6 (n)ll^ - ^ -1 ||i7(n) 



L 2 (n) 



We then integrate in time, apply several applications of the Holder inequal- 
ity, and utilize (5.16) to deduce 



<A*Eii<r 



1 10 

l£3(fi) 



0fc 



Uh — u 



k-l- 



L 2 {Q) 



+ CAt^2\\u k h \\\ 

k 

3 2 

< (At) 5 llefellloo^T;^^)) 



E A * 



2 

La(n). 



+Atc 



L6(Q) 



E K -Alli- 



en) 



< (At)lcrl + (At)^C ( ^ At 



2(7-3) 
5 7 



Ellefc-ejt- 1 ! 



5 7 



7 

L7(n) 



< C f (At)* + (At)*. 
This concludes the proof. 



□ 



Using the previous lemma, we are now ready to prove the following result 
and error bound. 



21 



KARPER 



Lemma 5.4. Let (Qh-> u h) be the numerical solution obtained through Defi- 
nition 3.1 and (3.7). Then, for all v £ L°°(0, T; VF 1,6 (f2)), 



T 

o Jn 



Dt(QhU h )I^v dx + ^2 I I \Jy>{qu ® u)U\v + dS(x)dt 



o J BE 



Dt(ghU h )v - g h u h <8> Uh ■ dxdt + 

where the reminder is bounded as 

i 

\F(v)\ < ^TC||Vu|| L oo( 0)T . L 3( n) ), 

unt/i constant independent of h and At. 

Proof. Using (5.2), we have the identity 
T 

o Jn 



(5.17) 



D^(g h u h )uXv dx + / \Jp(gu(g)u)Ulv + dS(x)dt 



o JdE 



f-i 4 />T 

/ Dt(g h u h )Il%v - g h u h (g>u h :X7v dxdt + V] / P («) d< 
Jo i=2 Jo 



Dt(QhUh)v - Q h u h ®u h :Vv dxdt 



4 /-T /.-j- 
+ V / Pi(v) dt+ D?(g h u h ) {lilv - v) dxdt. 

i=2 J ° J ° 



We thus define F(v) by 

4 „x 



\F{v)\ :- 



V / Pidt+ [ D^(g h u h )(ulv-v) dxdt 
„-_o Jo Jo 



i=2 
4 T 

Si 

4 r x 



r 


DtiQhUh) 


/o 





(SI) 



dt 



< V / |P*| (it + h(Aty 1 (At)^ \\Vv\\ L oo, 0T . 
i=2 J o 

4 r T i 
= V / \Pi\ dt + h^\\V 

;-i Jo 



£8(0)) 



v ll£°°(0,T;£3(Q))) 



where we have used Lemma 5.3 and h = a At. 

Finally, by applying (5.7) to bound the p terms, we obtain the desired 
(5.17). 

□ 



5.3. The artificial stabilization terms. To prove convergence of the nu- 
merical method we will need to prove that the artificial stabilization terms 
converges to zero. Moreover, we will need that these terms are small in a 
suitable Lebesgue space. 



A CONVERGENT METHOD FOR COMPRESSIBLE NAVIER-STOKES 



25 



Lemma 5.5. If(gh,Uh) is the numerical solution obtained by Definition 3.1 
and (3.7), the artificial stabilization terms satisfies 



^ e E 









r [ [Qhi T 


u Q h 4> 


dS(x)dt 


10 JT 







< h 12 C\\V(f>\\ L 2 i0tT . L 2 m , 

r-T 



JT 

13-6e 



Mr 



dS(x)dt 



< h 12 C[|Vt;|| £ oo(o ir; LS(n)), 
for all sufficiently smooth <p and v. 



(5.18) 



(5.19) 



Proof. We will begin by proving (5.18). By direct calculation, using the 
Cauchy-Schwartz inequality, we deduce 



^E 







[ L leklv 









dS(x)dt 



(5.20) 



= -/ i 1 - e V/ / Ia] r (n^) ds(x)dt 

= ~ hl ' e E / T / I^lr ( n h^ " A dS(x)dt 

J « c?-E7 

(E/ 7 / r Wr ^(x)dtj 
x (E f^h^-^f dS(x)dtj 

< h\h 1 ^ jf / r Wr dS(*)<ftJ l|V<^||i2 (0iT . i2(n))j 

where the last inequality comes from the trace theorem (Lemma 2.8) and 
the interpolation error estimate for LT^. Now, to bound the jump term, we 
apply Lemma 4.1 with B(g) = g 2 to obtain 

h l ~ e Y] [ T [ lQ h f r dS(x)dt 
p Jo Jt 



< I g dx + 

2 m 



T 



g\ div Uh dxdt 



io Jn 

^ lko||i2( n) + C\\g h 11^00(0^.^4(0)) || div Uh\\l^{0,T;L^{Q.)) 



(5.21) 



< 



1^0 1 1 L2( n ) 



8h\\L°°{0,T;L~i(Q.)) 



< ll£ ) o|li / 2(o) + h 6 C||£i/j||^oo( 0iT . L7 (o)) 
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where we have applied the inverse estimate (Lemma 2.9), the energy bound, 
and that 7 > 3. Finally, we apply (5.21) in (5.20) to discover 

r-T 



Mi 



Jr 



IT 



Q, 



dS(x)dt 



1 l-e 

< h^h 2 



^"E^/lrf dS(x)dt) 11^11^(0,2^(0)) 



< h 1 2 ^C\\V(f)\\ L 2 {0yT . L 2 {n)) , 
which is (5.18). 

Next, we prove (5.19). Since f r [n^f] r dS(x) = 0, we have the identity 

r-T 



r 




'0 Jr 

-T 



U- + u+ 
2 

U- + u 



Mi 



dS(x)dt 



JdE 



W r n[»-n[» ds(x). 



An application of the Holder inequality and several applications of the trace 
theorem (Lemma 2.8) allow us to deduce 

r-T 



Jr 



Mi 



dS(x)dt 



< 



h '~ e / r Wr dS ( x 

x sup ^ 



Iilv-U v h v 



T 



,)E 



U- + lt_|_ 



dS(x) 



dS(x) 



(5.22) 



^ hl ~ eC \Y,l / f Wr ^(*)j 

X (^^ll n h U - U h v \\L^(0,T;L3(n)) + ^^l|Vw|| iO o (0)T;L 3( n)) ^) 
X Ch~B ||% || L 2 (0,T;L e (n))- 

Next, we apply (5.21) and the fact that 7 > 3 to (5.22) to get 



Jr 



Mr 



dS(x)dt 



l-e ,21 1 



< /i 2 +3 6 iaC||Vu|| i oo( 0jT .jr ( S(n)) 
13-6e 

= h 12 C||Vv|| LO o( 0iT;i 3(n)), 



which is (5.19). 



□ 
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5.4. Weak time control. We conclude this section by proving /i-independent 
bounds on the discrete time derivates in the numerical method. 



Lemma 5.6. Under the conditions of the previous lemma, 

DtQ h E 6 Li(o,r ; w- 1 't(n)), 

D t {Q h u h ) £ h L\Q,T-W- l > 3 *{U)). 



(5.23) 
(5.24) 



Proof. Let <j) £ Cq°([0, T) x fl) be arbitrary and set as test-function in 
the continuity scheme (3.5) to obtain 



o ^si 



T 



Dt{Qh)<f> dxdt 



T 



Q, 



+ h 1 -* {g h j r 



n 



dS(x)dt 



o Jn 



(5.25) 



QhUhV(j) dxdt + P\(<j>) 



Whir 



n 



dS(x)dt, 



where the last equality is (5.1). From Proposition 5.2, have that 

|Pi(«)|<^q|v^|| 



'L 4 (0,T;Lir(r2)) 



(5.26) 



Moreover, from Lemma 5.5, we have that 

(■T 







[gfclr 



n 



r dS(x)dt < h^\\V4>\\ L2{0>T . L 2 m) . (5.27) 



Hence, by applying (5.26), (5.27), and the Holder inequality to (5.25), we 
conclude 



Dt{Qh)<t> dxdt 



o ^si 



< II^IIl2(O,T;L3(Q))||V0|| l2(OiT;L 3 (q)) 

+ /iiC7||V(/ ) || L4{0jT;L¥{Q)) + / l T||v ( A|| L 2 (0ir;L2(n)) . 

We conclude (5.23) by recalling that <fi was chosen arbitrarily. 

Next, let v G [Cq°([0, T) x Q)] d be an arbitrary vector and set Vh = n^f 
as test- function in the momentum scheme (3.6) to obtain 



T 

n Jn 

+ 



Dt(g h u h )v h dxdt + ^ 



o JdE 



Up (pit <8> u)vhdS(x)dt 



T 

o Jn 



VfcitftVfcVft - p(pfo) div^ u ft cfcccft 



d5(a;) = 0. 
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Now, to the first two terms we apply the identity in Lemma 5.4 and reorder 
terms to obtain 

rT 

Dt(QhUh)v dxdt 



o Jn 



/ / Qh^h ® Uh : Vv dxdt + F(v) 
Jo Jn 

+ / / ^hUhNhVh ~ p{Qh) div fe v h dxdt 
Jo Jn 

rT 



(5.28) 



Qhj dE Vh dS(x)dt. 



From Lemma 5.4, we also have the bound 

\F(v)\ < hic\\X7v\\ LoO ( 0iT . L 3 {n)) . 
Next, we invoke Lemma 5.5 to obtain the bound 

rT 



lQh\ dE Vh dS(x)dt 



(5.29) 



(5.30) 



< h 3 i2 6 C[| Vu||ioo(o t 



(0,T;L3(n))- 



By applying (5.29) and (5.30), together with the Holder inequality in (5.28), 
we deduce 



o Jn 



Dt(Qh u h)v dxdt 



< \\Qh\uh\ 



l L 2 (0,T;L2(Q)) 1 

+ h^C\\Vv\\ L oo { ^ T . L 3 {n)) 
+ l|p(efc)IU<-(o,TiLT(n))ll divv \\ L i {0tT . L i m 

13-6e 

+ h 12 C\\Vv\\ Lao{0jT . L 3 {n)) . 
From this we easily conclude (5.24). 



□ 



6. Higher integrability on the density 

In order to pass to the limit in the pressure, we will need higher integra- 
bility on the density. That is, from the energy estimate (Corollary 4.3) we 
only know that 

p(Qh) Gb L°°(0,T;L\n)), 

and L 1 is not weakly closed. Hence, it is not clear that p(gh) converges to 
an integrable function. To prove higher integrability, we shall make use of 
the operator A^-} : L p (ft) i-> W^ffl) 

A i [q]=-^A- l [U E q] , i = l,...,3. 

ax l q 

Here, He is the extension by zero operator to all of W 1 and |n denotes 
the restriction to Q. The A -1 operator is the usual convolution with the 
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Newtonian potential 

A-V = -A/ r*M-r dy, A>0. 
Jr3 \x - y\ 

Using Aii we define two operators A^ and .A dlv acting on scalars and vectors, 
respectively. For a scalar q and a vector v = [v±, V2, V3] T they are denned 

A S/ [q}= \A 2 [q] , A div [v]=A 1 [v 1 ]+A 2 [v 2 }+A 3 [v 3 ]. (6.1) 

By direct calculation, one easily verifies the following lemma. 

Lemma 6.1. For any two f G L p (n) and g G L«(0) mt/i 1 + i = 1 and 
1 < p, q < oo, the following identity holds 

[ vA v [g<p] dx = - I A dm [f\<Pg dx, G C °°(O). 

Moreover, there is a constant C such that, 

WAfihpw + l|v^[/]|| Z9(n) < c||/iu 9(n)) P < 

6 q 

We are now ready to prove higher integrability of the numerical density. 

Proposition 6.2. Let (gh,Uh) be the numerical approximation constructed 
through Definition 3.1 and (3.7). The following integrability estimate holds 

Proof. Let (f> G Cq°(Q) be arbitrary and define the test-functions 
V = <j)A V [Qh<l>]i v h = Il%v. 

Since Qh is piecewise constant in time, so is v and hence also Vh- Moreover, 
since <f> vanishes at the boundary, the degrees of freedom of Vh is zero at the 
boundary. As a consequence, Vh is a valid test-function in the momentum 
scheme (3.6). 

Note that the energy estimate, Lemma 6.1, and the Holder inequality 
provides the bound 

II Vu |U°°(0,T;L3>(fi)) < C||^||vFl,oc.(n)l|£'/ilU°°(0,T;L7(n)) < C||0|lw'i.«'(n)) ( 6 - 2 ) 

for any p < 7. 

Now, by applying Vh as test-function in (3.6), integrating in time, and 
reordering terms, we obtain 

/ / p{Qh) div/j v h dxdt = I 1 -\-I 2 + I 3 , (6.3) 
Jo Jn 

where 

h= I [ VhUhVhVh dxdt, 

iQhlr I^/ilr dS(x)dt, 



Jn 

T 
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'•:/' /• r'l /• 

Dt(QhU h )v h dxdt -V] / / Up(gu ® «) Ivfc] r dS(x)dt, 
io Jn r io ./r 



Before we start deriving bounds for Ii, I2, and ^3, let us first consider the 
term on the left-hand side of (6.3). For this purpose, recall from Section 2 
that the finite element spaces are chosen such that 

div^ Ii\v = njj divv. 

Hence, we have that 

p(g h ) div h v h dx = / p(g h )n® [div (4>A y [4>g h ])] dx 



p( 0h )dw(cl)A y ttQh}) dx (6.4) 
p(Qh)^4> • -A V [(t)Qh] + 4> 2 p{Qh)Qh dx. 
By setting this expression in (6.3), we see that 

/ / 4> 2 p{Qh)Qh dx = h+ h + h + h, 
Jo Jn 



(6.5) 



where now 



h = ~ / p(Qh)^4> ■ -4 V [#/i] dx. 
Jn 



Thus, the proof follows provided we can bound I\, I2, I3, and I4. 

1. Bounds on I\ and I2: Using the Cauchy- Schwartz inequality and the 
energy estimate (Corollary 4.3), we have that 



(6.6) 



\Il\ < ^^hUh\\L^(0,T,L 2 (n))\\^hVh\\L 2 (0,T-,L 2 (n)) 
< C\\Vv\\ L 2 {0:T . L 2 m < C||0||^ ljOo(n); 

where the last inequality is (6.2). 

To bound the I2 term, we apply Lemma 5.5 to obtain 

\h\ < /i5||Vu|| L oo (0iT;L7(c)) < fraC||0||^i,oo ( n)> ( 6 - 7 ) 

where again the last inequality is (6.2). 

2. Bound on I3: From Lemma 5.4, we have the identity 

h= f [ D?(g h u h )v - g h u h ®u h :Vv dxdt + F(v), (6.8) 
io Jn 

where F(v) is bounded by (5.17) and (6.2) as 

\F(v)\ < hlc\\Vv\\ L ^^ T . Lim < h^CU\\ 2 whoo{n y 
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It remains to bound the two other terms in I3. Let us begin with the easiest 
term. For this purpose, we apply the Holder inequality and (6.2) to deduce 
cT 

QhUh ® Uh : Vv dxdt 



Jn 

< \\QhUhU h \ 



L 2 (0,T;L~f- 1 (fl)) 1 



Wv\ 



|L°°(0,T;LT(Q)) (6-9) 

- i 1 1 ; 1 ' ~ 

L 2 (0,T;L^T(n)) 

where the last bound comes from Corollary 4.3 using that 



<CU\\ 



3 7 



> 



7 



for 7 > 3. 



3 + 7 7 - 1 

The remaining term in I3 is more complicated. By summation by parts in 
time followed by Lemma 6.1, we calculate 

cT 



Jn 



D t {QhUh)v dxdt 



M 



k=l 

M 



At 



v k dx 



k=l 



v k — v k 1 



At 



dx 



elulv 1 dx+ j Oj' f/,'' r Al r/.r. 

,k 



n 



To this identity, we apply the definition of v and Lemma 6.1 and write 

cT 

AI r 



Jn 



k=l 



M 



k=l 



A- ft 1 



g° h u° h v 1 dx+ I e ™u™v M dx 



At 

M~M„.M 



dx 



/ .4 div 




Jn 





At 



dx 



(6.10) 



glulv 1 dx+ I QhU^v M dx. 



M^MM 



The last two terms are easily bounded since (6.2) and 7 > 3 gives that 
v G b L°°(0,T; L°°(Q)) and hence 



Q° h vP h v l dx+ I Qfu%v M dx 



M~M„,M 



< C||0||| O o( n )||£i/ l 'U/ l ||i / o O (o ) r;ii(n))) 



where the last term is bounded by Corollary 4.3. To bound the other term 
in (6.10) we apply the continuity scheme (3.5) with 



qh = n. 



A 



div 
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M 



k=l 



/ .A div 




In 





o k - n k ~ l 
At 



M 



At EE [ u pV) 

k=i r 

M 

my: 



n 



Q 



A 



div 



d.r 



fc-i^fc-i 



dS(x) 



k=l 



g k u k V U div 



+ Pi [A 



i div 



Jc-l~fc-l 



da; 



where the last equality is (5.1) in Lemma 5.1. Now, by applying Proposition 
5.2 (i.e (5.7)) and the Holder inequality in the previous identity, we obtain 



M 



-AtJ2 / A 



k=l 



\ div 



o k - n k ~ l 
At 



dx 



< 2 ll^ tt fe||L 2 (0,T;L 2 (Q))ll</ > llwl, O o(Q) A dlV [Q h U h ] 

< 



L 2 {0,T;W^ 2 {n)) (6.11) 



2 



1 + hi ||^U/ l || L 4( 0; T ; L3(n)) 



where the last inequality follows from the properties of A" (Lemma 6.1) 
together with Corollary 4.3 giving gnUh Sb L 2 (0, T; L 2 (U)). Next, we apply 
the inverse estimate in Lemma 2.9 to bound the last term 

h*\\QhUhh*{o,T-,L3(n)) < ^ 2 (At)"4(7||£iA|| i 2( 0)T . i 3( n )), 

which is bounded by Corollary 4.3. By applying this in (6.11), using that 
h = a At, and setting the result into (6.10) we obtain 



T 



DtiQhUh)v dxdt 

lo Jn 

This, together with (6.9) and (6.8), yields 



<c\\<p\\ 2 wl ^ {n) . 



\h\ < 



(6.12) 



3. Bound on I4: To bound the I4, we will use that 7 > 3 yielding W 1,1 
embedded in L°° . The resulting calculation is 



\h\ = / p{Qh)^4> ■ A v [4>Qh] dx 
Jn 

< C||V^|| L °°(n)P V [#/ l ]||L-»(0,T;Vyi>7(Q)) < C\\V<p\\ 2 La0 ( Q) , 

where the last inequality is derived as in (6.2). Here, we have also used 
Corollary 4.3 which tell us that p(g h ) £ h L°°(0,T; L 1 (fi)). 
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4- Conclusion: By applying (6.6), (6.7), (6.12), and (6.13) in (6.5), we 
obtain 



'0 



P(Qh)Qh dxdt < CH^II^!,^^). 



Since p(gh) = agl and <p £ Co°(^) can be chosen arbitrary, this concludes 
the proof. 

□ 

7. Weak convergence 

In this section, we will pass to the limit in the numerical method and con- 
clude that the limit of {gh, u h) is almost a weak solution of the compressible 
Navier-Stokes equations. What remains in order to prove the main theorem 
is that p(gh) p{q)i which will be the topic of the ensuing sections. 

Our starting point is that Corollary 4.3 and Proposition 6.2 allow us to 
assert the existence of functions 

g € L°°(0, T; L 7 (f))), u e L 2 (0, T; W 1,2 (O)), 

and a subsequence hj — > such that 

Qh^g in L M (0,T;^(fi)), 

u h in L 2 (0,T;L 2 (O)), (7.1) 
V h u h ^Vu in L 2 (0,T;L 2 (fl)). 
From Corollary 4.3, we immediately obtain the integrability 
gu £ L°°(0,T;L m »(!]))nL 2 (0,r;L m2 (fi)), 
Qu®ue L°°(0, T; L\n)) n L 2 (0, T; L C2 (ft)), 

with exponents 

27 3 67 37 3 

^00 = — — > m 2 = — — > 3, c 2 = — — > -. 
7+12 3+7 3+72 

Moreover, Lemma 2.4 together with Lemma 5.6 provides the bounds 

g h e b C(0,T;Li(n)), Qh u h G b C(0,T;L m ™(Q)). 

In the following lemma, we prove that the above convergences are suffi- 
cient to pass to the limit in both nonlinear terms involving u. 

Lemma 7.1. Given the convergences (7.1), 

QhUh, QhUh QU in L 2 (0, T; L m2 (JT)), 

g h u h - qu in C(0,T;L m °°(ft)). (7.2) 

QhUh ®Uh—^QU®uin L 2 (0, T; L C2 (O)), 

Proof. Lemma 2.10 tell us that Uh is spatially compact in L 2 (0, T; L P (S7)) 
for any p < 6. From Lemma 5.6 we have that DiQh Gb -^ 5 (0, T; W" 1 '^ (f2)). 
As a consequence, we can apply Lemma 2.6, with = gh and fh = Uh, to 
obtain 

g h u h ^gu in L 2 (0, T; L m2 (Q,)). 
To conclude weak convergence of QhUh and QhUh, we write 

ghUh = ghUh + Qh(P-h u h - u h), Qh^h = ghUh + ghij^u h - u h ). 
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The interpolation error estimate (Lemma 2.7) tell us that 

\\^ u h-Uh\\L 2 (O,T-L 2 (a)) + \\^-h U h~ u h\\L 2 (0,T;L 2 {Q)) < C\\^ h,Uh\\L 2 (0,T;L 2 {n)), 

and hence QhUh — 1 f?/i«/i — 1 in the sense of distributions on (0, T) x Q. 

Next, since D^(g h u h ) <G b L 1 (0,T;T^- 1 '2(O)) (from Lemma 5.6), another 
application of Lemma 2.6, this time with = QhUh and fh = u k> renders 

QhUh®u h -± gu®u in L 2 (0,T;L C2 (ft)). 

Clearly, this also implies that Qh u h &> % — ^ qu (3 u. The convergence of 
f?ft^ft Sfe then follows by writing 

Q h U h ®U h = QhUh ®Uh + Qh{^h U h ~ Uh) ® Ufc, 

and applying the interpolation error estimate on the remainder. 

□ 

Using the convergences we have derived thus far, we are able to pass to 
the limit in the continuity approximation (3.5). 

Lemma 7.2. Let (gh,Uh) be the numerical solution obtained by Definition 
3.1 and (3.7). The limit (g,u) is a weak solution of continuity equation: 

Qt + div(gu) = in V'([0,T) x H). 

Proof. Let <$> G Cq°([0, T) x ft) be arbitrary and set as test-function in 
the continuity scheme (3.5) to obtain 
i-T r 

Dt(g h )4> dxdt 



o Jn 



E 



T 







n 



Q, 



T 

o Jn 



XJp(gu) 
Qh,UhS74> dxdt + P\{4>) 



where the last equality is Lemma 5.1. 
From Proposition 5.2, we have that 

i 



n; 



dS(x)dt 



(7.3) 







fo / r Wr 









dS(x)dt, 



\Pi((f>)\<h*C\\cj>\\ 



L 4 (0,T;LTT(n))' 



(7.4) 



From Lemma 5.5, we have the bound 

T 



E 



o 



^ £ Wr 



IT 



dS(x)dt 



<h-^C\\Vct>\\ L 2^ T . L 2 m . (7.5) 



Summation by parts provides the identity 

f [ D^Qh)^ dxdt 

Jo Jn 

rT—At r 







n 



g h Dt((fi(+At)) dxdt - / Q° h <f)(At) dx 



(7.6) 
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Now, we apply (7.6) in (7.3) and send h — >• along the subsequence where 
QhUh — Qu and apply (7.4) and (7.5) to obtain 

g(4>t + u • V0) dxdt = — I £o</>(0, •) dx. 



□ 



/o Jn 

This concludes the proof. 



Next, we prove that the limit (g, u) is almost a weak solution of the 
momentum equation (1.2) and hence that it only remains to prove strong 
convergence of the density. 

Lemma 7.3. Let (gh,Uh) be as in the previous lemma. The limit (g,u) 
satisfies 



(gu) t + div(gu ® u) + Vp(g) - Au = in T>'([0,T) x Q), (7.7) 



where p(g) is the weak limit of p(gh)- 

Proof. Let v € [C^°([0, T) x $7)] 3 be arbitrary and set Vh = U^v as test- 
function in the momentum scheme (3.6). After an application of Lemma 
5.4, we obtain the identity 

T 

Dt(QhUh)v - ghUh Uh ■ dxdt + F(v) 



o JCl 

+ / / V h u h V h Vh - p{gh) div h v h dxdt (7.8) 
Jo JQ 

Lemma 5.4 also provides the bound 

\F(v)\ < fcM|Vu|U» (0)T . £ 3 (n)) . (7-9) 
Moreover, from Lemma 5.5, we have the bound 



r 

< /l^C||Vu|| L oo (0iT;L 3 (n)) . 



(7.10) 



By applying summation by parts to the first term in (7.8), sending h — > 
along the subsequence for which ghUh ® — v gu ® u (Lemma 7.1), and 
applying (7.9) and (7.10) we obtain 



/ / —guvt — gu®u : Vv + VuVt> — div v dxdt = 
Jo Jn Jo, 

which concludes our proof. 



mov(0, •) dx, 



□ 
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8. The discrete Laplace operator 



In the next section we establish the most important ingredient in our 
proof of compactness of the density approximation, namely weak continuity 
of the effective viscous flux. However, before we can embark on this proof, 
we will need to establish some properties related to our numerical Laplace 
operator. 

In contrast to a standard continuous approximation scheme, our discrete 
Laplace operator does not respect Hodge decompositions. More specifically, 
to prove the upcoming Proposition 9.1 we shall need to use test-functions 
of the form 



where we for the purpose of this discussion does not require v to satisfy 
boundary conditions. At the continuous level, testing with this v is equiva- 
lent to applying „4 dlv [-] to the equation. In particular, since curl„4 v = 0, v 
satisfies 



This property is not true for our discrete Laplace operator. However, in the 
upcoming analysis it is essential that, at least, 



This is the result that we will prove in this section. 

As will become evident, the property (8.1) is not trivially satisfied for our 
discretization. It is the extra "artifical" stabilization terms in the scheme 
(adding diffusion in all directions) that will provide the needed ingredient. 
In fact, the property (8.1) is the only reason for the presence of these terms. 
This discretization strategy was devised by Eymard et. al for the stationary 
compressible Stokes equations [9]. In fact, most of the material contained 
in this section can be found there with only slight modifications to fit the 
present case. 

In our proof, we shall need the operator : Qh{&) Pft(^) interpolat- 
ing piecewise constant functions in the space of continuous piecewise linears 
Vh (Lagrange element space). This operator is defined by 



for all vertices r in the discretization, where N T is the set of elements having 
r as a vertex. Note that shape regularity of Eh renders the cardinality of 
N T bounded. The following result can be found in [9, Lemma 5.8]: 

Lemma 8.1. Let qh £ Q^(fi). There exists a constant C , depending only 
on the shape-regularity of Eh such that 



[Q] 





a > 0. (8.1) 





(8.2) 




(8.3) 
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We shall need the following auxiliary result. 
Lemma 8.2. Let Uh G V^(O) and v £ T4 /2,2 (0) &e arbitrary. Then, 

/ VhUh^vdx= / curl/j n/j curl/j t; + div^ Uh divv dx + E(v, Uh). (8.4) 
Jn Jn 

Furthermore, there is a constant C > 0, depending only on the shape- 
regularity of Eh such that 

\E(v,u h )\ < hC\\V h u h \\ L 2( Q) \\X7 2 v\\ L 2( n) . 

Proof. Using the Stoke's theorem and the identity —A = curl curl — V div, 

= 2. / —UhAv dxdt + / (Vi> • v)uh dS{ 
p Je Jo JdE 



x ) 

T 



y / Uh curl curl v — UhV div v dxdt + / / 
E Je Jo Jd 

curl/j Uh curl v + div Uh div v dx 



(Vv ■ v)uh dS(x) 

dE 



n 



+ 2_. / (curly x u)uh x v — (div v)uh ■ v + (X7v ■ v)uh dS(x), 
E JdE 

which is (8.4) with 

E(v) = / (curly x v)uh x v — (divu)u^ • v + (Vv ■ v)uh dS(x). 

g JdE 

Next, we use that the trace of Vv is continuous across edges to deduce 
\E{v)\ < V I \Vv\u h dS(x) 

E JdE 

= J2 [ \Vv\(uh-ltfu h ) dS(x) 

= Y, I l Vu - U h (V«)|(« A - U^uh) dS(x) 

^ JdE 



E 

1 



< h-ic [\\vv - n«(W)|U 2(n) + h\\v 2 v\\ L 2 (n) 

x ^ (\\ Uh - n®u h \\ L 2( n) + ^HVfeU^H^^)^ 

< hC\\V 2 v\\ L 2 {n) \\\\Vu h \\ L 2 (n) , 

where we have applied the trace theorem (Lemma 2.8) and the interpolation 
error estimate (Lemma 2.7). This concludes the proof. □ 

As in [9], we obtain the following identity. 
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Lemma 8.3. Let (gh,Uh) be the numerical solution obtained by Definition 
3.1 and let A v [) be given by (6.1). For any ip G C°°([0,T) x Q) and <f> e 



[ V h u h V h I$ (iPA V [<t>Qh}) dxdt 
o Jn 

= / curl fe u h X7ip x A V [4>Qh\ + div h u h ■ A V [4>Q h ] dxdt 
Jo Jn 

+ / / <pip drvhUfiQh dxdt + F{gh), 
Jo Jn 



where the T term is bounded as 

\HQh)\<h^-^C,, (8.5) 

with C independent of h and where we recall the requirement e > g . 
Proof. To simplify notation, let 

v = ^ [<p eh ] , v L = [<t>n% Qh ] , 

and observe that linearity of A^ provides the identity 

v - v L = ^A v [(p{Q h -U^g h )). 
By using Lemma 2.11 and adding and subtracting vl, we deduce 



r 

o Jn 



VfeWftVfell^ dxdt 



T f 

/ V/jM^Vf dxdt 

o 



/ VhUh^VL dxdt + / VhUhV(v - vl) dxdt 
Jn Jo Jn 

T f - 

curl/j u/j curl + div^ Uh divvL dxdt 



o 



o Jn 

T 



+ 



/ / V/ l u /l V(u - vl) dxdt + E(v L ,u h ) 
Jo Jn 
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where the last equality is Lemma 8.2. Next, we once more add and subtract 
vi to obtain 
cT r 

^hUh^h^iv dxdt 

o Jn 

T f 

/ curl/j Uh curl v + div/j Uh div v dxdt 
ii Jn 

T 



+ 



+ 
T 

u Jn 

rT 



/ / curl^ Uh curl(i;£, — v) + div^ Uh &\v{vl — v) dxdt 
Jo Jn 

/ / V h u h V(v - v L ) dxdt + E(v L ,u h ) 
Jo Jn 

curl ft u h Vifj x A v [4>Qh] + div/,, UhVip ■ A v [4>Qh] dxdt 



1.6) 



+ / / 4>ifjdw h u h g h dxdt + ^(qh), 
Jo Jn 



where we have used the definition of vi and introduced the quantity 

F{Qh) ■= / curl^ u h eurl(i>£ — v) + div^ uu d\v{vL - v) dxdt 
Jo Jn 

\ \ V h u h V(v - v L ) dxdt + E(v L ,u h ). 
Jo Jn 

In view of (8.6), it only remains to prove (8.5). 

Some applications of the Cauchy-Schwartz inequality together with the 
properties of A^ 7 [■] gives 

L 2 (0,T;L 2 (B)) + \E(v L ,U h )\ 

l 

-T c IT T 2 \ 2 



hC \Y,f [^dS(x)dt) +h\\V 2 v\\L*(0,T;mn)) 

\ TcB Jo Jr J 

< Vh2C t V [ T [ lg h j z r dS(x)dt 
\ rcB Jo Jr 



where the second last inequality is (8.3) and Lemma 8.2 and the last in- 
equality is and application of (8.3). To bound the last term, we invoke 
(5.21) yielding 

\J"(g h )\ < h^-^C, 

which is (8.5) 

□ 



9. The effective viscous flux 

The main tool that will allow us to conclude strong convergence of the 
density is a remarkable result discovered by P. L. Lions [ ] for a continuous 
approximation scheme. The result states that the effective viscous flux 



divu — p(g), 
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behaves as if it is converging strongly. More specific, in our numerical set- 
ting, the result goes as follows: 

Proposition 9.1. Let (gh, Uh) be the numerical solution obtained using Def- 
inition 3.1 and (3.7). Moreover, let (g,u) be given through the convergences 
(7.1). Then, 



lim / / 4>ifj (div h u h - p(g h )) g h dxdt 
Jo Jn 

= J J rfnp (div u — p(g) S J g dxdt, 

for all <f> £ Cg°(fi) and ip £ C§°([0,T) x fi). 

Hence, the product of the weakly converging effective viscous flux and the 
weakly converging density converges to the product of the weak limits. Our 
proof of this proposition will rely on a number of auxiliary results which we 
will prove first. The proof is concluded in Section 9.2. 

9.1. The numerical commutator estimate. In the upcoming analysis, 
the following lemma will be essential. A proof based on the div-curl lemma 
can be found in [10]. 

Lemma 9.2. Let v n and w n be sequences of vector valued functions such that 
v n — v in L P {Q.) and w n — 1 w in L q (Q), respectively, where 1 < p, q < oo 
and i + | < 1. Moreover, let B n — B in L p (£l). Then, 

(1) v n VA div [w n ] - w n VA dw [v n ] — >• vVA dw [w]-wVA dw [v] 

(2) B n X7A dm [w n ] -w n VA v [B n ] — > BVA dm [w] - wVA v [B] 
in the sense of distributions on fL 

Lemma 9.3. Given the convergences (7.1) - (7.2), 



'0 Jn 

Proof. We begin by observing the following identity 

l-T 



lim / / <pg h u h V (A dm [ipg h u h ] ) - ipg h Uh ® u h : V (A v [4>g h ]) dxdt 
(pguV (A dm [^gu] \ -i/;gu®u:V {A^[(/>g}) dxdt. 



o Jn 

T 



o Jn 

fT 



4>ghU h V \ A dlv [tpQhUh]) - ipghUh 0u h :V (A^[4>g h ]) dxdt 
J' J^h- [V (A d ™[^g h u h \) 4>Qh ~ V (A v [cf>g h ]) ■ ^g h u h ] dxdt (9.1) 



UhH h dxdt, 

o Jn 

h 



where we have introduced T~L given by 

H h = V U div l?Pg h u h ]) <f>Q h - V {A*[<t>Q h ]) ■ %bg h u h . 
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From Lemma 7.1, we have that 

g h ^ gin C(0, T; L 7 (fi)) n L 2 (0, TL 7 (ft)), 

27 (9.2) 
QhU h in C(0, T; (0)) n L 2 (0, T; L 3 (0)). 

Since in addition 7 > 3, the Holder inequality, Lemma 6.1, and Corollary 
4.3, provides the estimate 

11% 11^2(0 T-ii(n)) - c 'll£ >h ^lli 2 (o,T , ;Z,3(s^)) || £>^|Uo°(o,r;iT(^)) < C, (9.3) 
and similarly, 

11% II , < C||^W/i|| , -2r, Jl0/JU°°(O,T;.t/r(n)) < C, 

(9.4) 

where > 1 since 7 > 3. 

By virtue of (9.2)-(9.4), we can apply Lemma 9.2, with vh = ipQhUh and 
Bh = 4>Qhi to conclude the weak convergence 



U h = V U dlv [^£>/A] <t>Q h - V (-4 V [MJ) • ^Q h u h 

/ \ (9-5) 

V (A dlv [tpQu}) (pg-V (A v [<t)Q\) ■ ipgu =: H, 



in L 2 (0, T; L§ (O)) n C(0, T; L^< (SI)) as ft -> 0. 

To proceed we will need the standard mollifier which we will denote by 
i?* 5 . It will be a standing assumption throughout that 5 is sufficiently small. 

Now, by adding and subtracting we write 

rT f 

/ u h H h dxdt 

^ rp (9.6) 

= f f (u h - u h )U h + {u h -R 5 * u h )H h + (R 5 * u h )U h dxdt 
Jo Jn 

The first term in (9.6) converges to zero as 

[ f (uh- u h )U h dxdt 
Jo Jn 

< \\uh ~ u h \\ L 2 { ^ T . L3m WH h \\ 3 



< /i||V/ l u/ l || L 2( 0iT;L 3(n))||^ h || 2 3 



L2(0,T;L2(n)) (9 7 ) 
L 2 (0,T;L2(C)) 

< ^||V,n,|| L2{0iT;L2(n)) ||^|| i2(0)T;L 3 (Q)) < hIC. 

To bound the second term in (9.6), we apply the Holder inequality and the 
space translation estimate of Lemma 2.10, 

fT 



f [ (u h -R 5 * u h )U h dxdt 
Jn 



< 



u h - R s * u h 



\\U h \\ 3 (9-8) 



< C {h 2 + 5 2 ) 3 || V h U h \\ L 2 (0jT . L 2 m <C{h 2 + 5 2 ) 3 . 
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■2; 



Next, since U h in C(0,T;L^(ft)), we have in particular that 



% in L 2 (0,T;W~ 1,p (tt)), p< 



Hence, for each fixed 5, we can conclude that 



lim 



o Jn 



(R s *u h )H h dxdt 



o Jn 



(R s * u)U dxdt. 



Finally, we pass to the limit in (9.6) and apply (9.7) and (9.8) to obtain 



lim 

7i->0 



o Jn 



UhH h dxdt 



o Jn 



(R s *u)n dxdt + 0(V5). 



We conclude the proof by sending (5—^0 and recalling (9.1) and (9.5). 



□ 



Lemma 9.4. Let (gh,Uh) be the numerical solution obtained through Def- 
inition 3.1 and (3.7). Let (g,u) be the corresponding weak limit pair ob- 
tained from the convergences (7.1). Then, for any <p G Cq°(£1) and ip S 
C o °°([0,T)xfi), 



lim > 

/w0 > ^ 



r 

o Jr 



o Jn 



-\Jp(gu®u) itflll [^ V [#h]]J r dS{x)dt 
guV{c^A dw {^gu)) - gu®u:V {^A V [4>g]) dxdt, 



where the operators A^[-] and A dw [-] are given by (6.1). 

Proof. Our starting point is provided by Lemma 5.1 which in this case gives 



o Jr 



Up(fm) [^ div [V^] 

- \Jp(gu u) lu^Ul [tpA v [(j>Q h }) dS(x)dt 



o Jn 



ghUh^ [4>A lv [ipghUh]J - QhUh ® % : V (i/>.4 dxdt 

4 

+ Pi (<M div [Vw]) + ^p (M v [&?/.]]) • 



(9.9) 



i=2 

Let us first prove that the Pi terms are converging to zero as h — > 0. 
1. Convergence of Pi to zero: From (5.3), we have the following bound 

P 1 {<t)A^[iPg h u h }) 

' V (cjiA^&QhUh]) (9.10) 



1 ,4- 7 
< h* ^ 47 C 



L 4 (0,T;L3(Q)) 



1 9H 



Z, 4 (0,T;W 1 - 3 («)) 
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From the regularization property of ^4 dlv (Lemma 6.1) together with an 
inverse estimate in time (Lemma 2.9) we deduce 

\A diY [ipQ h u h ] 

< C\\tpQ h U h \\ L A^ !T . L 3^ , g ^ 

< C||^|| L oo(( 0iT ) xQ )||^/ l n/j|| L 4( 0jT;L 3(Q)) 

< (A) _ 4C||^|| LOO((0jT)xC) || £ ) /l n/j|| i 2( 0>T;i 3(n)), 

where the last term is bounded by the energy (Corollary 4.3). By setting 
(9.11) in (9.10) and using that At = ah, we obtain 

P l (^A^^QhUh]) I < h*- 3 ^C, (9.12) 

where the exponent for h is strictly positive as 7 > 3. 

For the remaining P l terms, we apply (5.7) in Lemma 5.1 to obtain 

4 

j2\Pi(M v [^h}})\ 

(9.13) 

< /l 3 C||V'||Loo( 0)T . H /l,oo( n ))||(/)|| i cx)(Q)||£)/ l || LO o( 0jr . i7 ( n )) < hiC, 

where we in the second last inequality have applied Lemma 6.1. The last 
inequality follows from p(g h ) G b L°°(0,T; L 1 (O)). 

Consequently, (9.12) and (9.13) allow us to conclude that 

4 \ 



i=2 



<h\c\\V{iPA W [<pQh]])\\ 



lim 



A UA div [i>QhU h })\ +Y,\ P i (M V [#/,]]) I = 0. (9.14) 



i=2 



2. Convergence of the commutator term: To prove convergence of the first 
term after the equality in (9.9) we will need to rewrite this term on a form 
for which Lemma 9.3 is applicable: 

g h u h V (4>A dlv [ipQ h u h ] ) - g h u h ®u h :V (ipA v [(j)Qh\) dxalt 



T 



Jn 

VUh ah v [ ^ 

Jn 



(j)QhUhS7 [A lv [i/jQhUh]J - i^QhUh ®u h :V [A [(j>Qh\) dxalt 

(9.15) 

+ / / QhUhVc/) ■ A A "[i>QhUh] - QhUh ®u h :V'4)(3 A V [4>Qh] dxdt. 
Jo Jn 

Now, observe that the first term after the equality is precisely the one covered 
by Lemma 9.3. Hence, 

lim / / 4>g h u h V \A dlv [^Q h u h ] ) - i^QhUh <8> u h : V (A v [4>g h ]) dxdt 

ft->0 /n let V / 



Jn 

T 



</)quV (A div [4)Qu] ) - ifjgu ®u:V (A^ 7 [4>q]) dxdt. 
Jn v ' 

(9.16) 
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Moreover, in view of the convergences in Lemma 7.1, we have that 
A v [<t)Q h ] -> A v [^q] in (7(0, T; L P {£1)) for any p < oo 
A^Q h u h ] -)• .4 div (V^u) in C(0, T; L 3 (ft)). 
As a consequence, there is no problems with concluding that 

lim / / g h u h X7(f> ■ A diy [ipQ h u h ] - g h u h (g> u h : <g> A v [(/>Q h ] dxdt 
h^oJo Jn 

T r 

/ guV(f) ■ A d ™{4>Qu] - Qu®u:Vif)®A v [<f>o\ dxdt. 
o Jn 

(9.17) 

Finally, we send h — > in (9.15) and apply (9.16) and (9.17) to conclude 
lim / / Qh,u h V \ $A diy [ipQhUh]) ~ QhUh ®u h :V (ipA^ '[<j)Qh]) dxdt 



r 

o Jn 



cf)QuV [A dlv [tpgu]J - ipgu 0u:V (A^[4>q]) dxdt 

T 



guV(f> ■ A aiv [i[>Qu] - gu®u:Vip®A w [4>q] dxdt (9.18) 



T 

o Jn 



div r 

yu v yj • ^-i 

(I Jn 

guV{il)A d "{4>Qu)) - qu®u:V (tpA V [(j)Q}) dxdt. 



3. Conclusion: We conclude the proof by sending h — > in (9.9) and 
applying (9.14) and (9.18). 

□ 

9.2. Proof of Proposition 9.1. Define the test-function 
v = ipA V [(f>Qh], v h = Il%v. 

By setting as test-function in the momentum scheme (3.6), integrating in 
time, applying Lemma 8.3 for the term involving V 'hUh\ 'hVh, the calculation 
(6.4) for the term involving the pressure, and reordering terms, we obtain 

/ / H{p(Qh) ~ div h u h )Q h dxdt = Vjf + F(Qh), (9-19) 
jo Jn ,. =1 

where F(Qh) is given by (8.5) and 

J i=~ f f P(Qh)V<t> • A v [(t>Q h ] dxdt, 
Jo Jn 

J$= / curlfc u h Vip x A v [<j>Q h ] + div h u h Vip ■ A v [4>Q h ] dxdt, 
Jo Jn 

4=1 I D^{g h u h )v h dxdt-Y] [ [ Up(eu ® S) [u & ] r dS(x)dt, 
Jo Jn r Jo Jt 



r 

Jh = h 1 -* £ J v {^Y^) IMr Mr dS{x)dt 
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Next, define the test-function 

w = ipA v [M- 

Observe that this test-function is the limit of v. That is, 

v — > w a.e as h — > 0. 

Setting w as test-function in the weak limit of the momentum scheme 
(7.7), and reordering terms, yields 

f ^(p(pJs) - divu )Q dxdt = Ji + J 2 + J 3 . (9-20) 

T 



Ji = - / p(g)V(f>- A v [<fig\ dxdt, 
Jo Jn 

J 2 = / curl uVip x A v [4>q\ + div uVip ■ A v {(fig] dxdt 
Jo Jn 

J3 = — / / guwt + gu®u : Vw dxdt — / mow(0, •) dx. 
Jo Jn Jn 

Now, observe that the proof of Proposition 9.1 is complete once we prove 

Since we have already established the convergences 

-4 V [#ft] -»■ -4 V [#] in C(0, T; L p (ft)) for any p < 00 
^[^fc] -»• -4 div (V>H in C(0, T; L 3 (ft)), 
we can immediately conclude that 

-> Ji J$ -»■ J 2 . (9.21) 
For the J4 term, the calculation (6.7) gives 

|J 4 h | < h*C->0. (9.22) 

Hence, the only remaining ingredient in the proof of Proposition 9.1, is to 
establish J3 — > J3. Indeed, let us for the moment take this result as granted 
(Lemma 9.5 below). Then, we can send h — > in (9.19), using (8.5), (9.21), 
(9.22), to discover 

lim / / 4)il)(p(g h ) - A\\ h u h )g h dxdt 
Jo Jn 

4 

= lim + F{g h ) = Ji + J 2 + h 



i=l 



Jn 



4>ip(p(g) — div/j u)g dxdt, 



where the last equality is (9.20). This is precisely Proposition 9.1. 
Hence, the proof is complete once we establish the following lemma. 
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Lemma 9.5. Under the conditions of Proposition 9.1, 



lim4 = J 3 , 

where J3 and J3 are given by (9.19) and (9.20), respectively. 

Proof of Lemma 9.5. To prove convergence of J3 , we shall need to rewrite 
the time derivative term in J3 1 using the continuity scheme. By adding and 
subtracting, and applying summation by parts, we deduce 

r-T r 

Dt(g h u h )v h dxdt 



Jn 



T—At 



g h u h D?(v(+At)) dxdt 

■T 

,0^0 



QtuMAt) dx + 



T 

At Jn 



Jn 

{g h u h ){- - At) D h t (^ v [<f>g h ]) dxdt 



Dt(g h u h )(v h - v) dxdt 



g° h u° h v(At) dx + 



At Jn 



( 6h u h i>)(- - At) A^ 



Jn 



Dt(g h u h )(v h - v) dxdt 



dxdt 



T 



(g h u h )(- - At) A v [</>g h ] D t V dxdt 
T 



Q%Sy h v(At) dx + 



Jn 



Dt{g h u h )(v h - v) dxdt 



L>l,iil,)\- - _k ) , jX ' - 1 r ' ,h ■ 

1 At Jn 
in 

Next, we apply the integration by parts formula for A^[-] (Lemma 6.1) 
i-t r 

'o Jn 



D't (QhUh)v h dxdt 

7 

'Ai Jn 



A div [(Q h u h iP)(--At)]<pD? 6h 



[ [ (Q h u h )(- - At) A v [</>g h ] dxdt 
J At Jn 



Q u h utv(At) dx + 



T 



Jn 



Dt(g h u h )(v h - v) dxdt. 



We then rewrite the first term using the continuity scheme (3.5) with 



qh = n,_ 



ct>A^[{g h u h ^)(- - At)} 



The resulting expression reads 
(■T 

Dt(g h u h )v h dxdt 



Jn 



At Jr 



Up(gu) 



ir 



Q 



(t)A A "\{g h u h ip){- - At)) dS(x)dt 
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cj>A div [(g h u h ^)(- - At)} dS{x)dt 



At Jfl 
0^ 



{g h u h ){- - At) A v [4>Q h ] Dfy dxdt 

T 



(9.23) 



Qlulv(At) dx + / D? (QhUh)(vh - v) dxdt, 
n Jo Jn 



By applying (9.23) in the expression for Jg , we obtain 

4 = E [ T [ Up(eu) fn? [M div [(^v)(- - At)]l] ^(x)dt 

p iAf Jr JJIr 

- E I I_Vv{qu®u) fnjitf [g h <f>]] 1 ^ dS(x)rfi 



o Jr 



^ £ E 



o 



[eh] 



h 



<t>A d "[{g h u h ^)(- - At)} dS(x)dt 



At JO, 
0^0 



(g h u h )(- - At) A v [4>Q h ] D?4> dxdt 

T 



g u h u v h v(At) dx + I I D?{g h u h )(v h -v) dxdt. 
n Jo Jo. 



In order to apply Lemma 9.4 we need to remove the time shift —At in the 
first term. By adding and subtracting, we obtain 



4=E 



Vp(gu) 11® cj)A d "[QhU h if>] dS(x)dt 
At Jr II JJIr 



\Jp(gu ® u) fnjjnjf [ipA^ [Qh<t>]] 1 dS(x)dt 



o Jr 

T 



E 

r 

E 

r 

^ £ E 



Up(e«)At n£ (PA d "[Dt(ghU h ^)} dS{x)dt 
At Jr 11 Jllr 



At Jr 



[Qhl 



n 



h 



4>A div [(QhuM- - At)] dS(x)dt 



T 



At Jtt 
0^0 



(g h u h )(- - At) A v [<t>Q h ] Dty dxdt 

T 



(9.24) 



Qlv? h v(At) dx+ / D?(g h u h ){v h - v) dxdt. 
n Jo Jn 



i=i 



where the third term contains the time difference. 
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Now, observe that K\ + K 2 is exactly the terms covered by Lemma 9.4. 
In particular, 



lim K\ + K 2 



JS1 



guV 4>A dlv [ipQu] - gu®u: V(^A V [<pg]) dxdt. 



(9.25) 



To bound the K3 term (the new term), we apply Lemma 5.1 and Propo- 
sition 5.2 to deduce 



\K» 



E 



Jr 



Up(gu)At ct)A div [D%(g h u h ip)] dS{x)dt 



At 



T 



At Jn 



g h u h : V UA div [D?{Q h u h il>)] ) dxdt 



+ (At)P 1 ( K M div lDHQhU h iP)} 
67 \\ghUh(t) - Qh,v 

^' 0} C\\g h u h (t) - g h u h (t - At)\\ 



(9.26) 



<\\QhUh\\ 67 \\g h u h (t) - g h u h (t - At)\\ e 7 

L 2 (0,T;L3+^(C)) L 2 (At,T;L^f^(n)) 



1 max{3 i^, } 



+ h 2 

< C\\g h u h (t) - g h u h {t - At)\\ 



L 4 (At,T;L~5~ (Q)) 



67 



L 2 (At,T;L5^(tr)) 

+ ^C|| efc « fc (t) - eh « fc (f - A*)ll L2 _ a(At)T;i¥(n)) 

where 0(7, 5) = | — max{3^^,0} — 2 (2-&) anc ^ <5 > is a small number. 
To achieve the previous inequality, we have applied the inverse inequality in 
time (Lemma 2.9). 

Now, from Lemma 5.3, we know what 

QhMh(t) — QhUh(t — At) — > a.e on (0, T) x Q. 

Hence, since g h u h e L°°(0, T; L§ (0)) n L 2 (0, T; L 3 (ft)), 

ff/,«fc(t) " Q h u h {t-At) ^ in L M (0.r;i ?1 (^)nI P2 (0,r;L« 2 (fi)), 



for any p\ < 00, gi < , P2 < 2, and q 2 < jpn. Moreover, since 7 > 3, 
can choose 5 small such that 0(7, 5) > 0. As a consequence, we can pass to 
the limit in (9.26) to obtain 



67 



wc 



lim \K 3 \ = 0. 
The K4 term is directly bounded by Lemma 5.5; 

i-T 



(9.27) 



\Ka\ 



At Jr 



n 



h 



c/)A dW [g h u h ij)(- - At)] dS(x)dt 



< h^C I V (<f>A dW [QhU h i>{- - At)} 

ll-4e ^ 

< h s C 1 1 gfetih 1 1 x,a (o jT;i 2 ( n , 
where the norm is bounded by Corollary 4.3. 



L 2 (At,T;L 2 (n)) 



(9.28) 
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Next, since 

Q h u h ^QU inL 2 (0,T;L m2 (ft)),m 2 >3, 
and Ji 7 [<j)Qh] — > A^[4>g] i n C(0, T; L p (fl)), for any p < oo, we see that 
lim (K 5 + K 6 ) 

h^-0 



lim / / (g h u h )(- - At) A v [<t>Q h \ Dftp dxdt 
J At J n 

r (9 29) 

lim / g° h u? h v(At) dx 



guA^ [4>g]ipt dxdt — / mQWQ dx. 
/o Jn Jn 

To bound Kf, we apply Lemma 5.3 and the interpolation error for LT^ 
and obtain 

\K 7 \= / / Dt(g h u h )(v h - v) dxdt < hi C\\Vv\\ L aor 0T . L 3r n )) , 

Jo Jn (9.30) 

JL 

< hiC\\gh\\L°°(o,T;Li(n)), 

which is bounded by Corollary 4.3. 

At this stage, we can send h in (9.24), using (9.25), (9.27), (9.28), 
(9.29), and (9.30), to obtain 

lim 4=1 [ guV (<j)A div [ipQu]) - gu®u: V(^A v [(pg}) dxdt 
Jo Jn v ' 



T 

o Jn 



guA^ '[(f>g]ipt dxdt — / rriQWo dx. (9.31) 



Finally, since the limit (g, u) is a solution to the continuity equation 
(Lemma 7.2), we have that (recall that c/> does not depend on time) 



l T f guV UA d "[^gu] 
'o Jn v 



<j)g(A div [ipgu]) dxdt - / 0g o A div [ip(O, -)m ] dx 
o Jn '* Jn 

T 

(A^[g(J)]) t tl>gu dxdt 

guwt dxdt + I I gu A^ 7 [4>g] ifjt dxdt. 
Jo Jn 



o Jn 

T r 



Jo Jn 

Thus, by using this identity in (9.31), we obtain 

lim J3 = — / / guwt + gu® u : V«i dxdt — I mow(0) dx 



h^O 

= J 3 ; 



Jn 



which was the desired result. □ 
The proof of Proposition 9.1 is now complete. 
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10. Strong convergence and proof of Theorem 3.5 

Equipped with Proposition 9.1 we can now obtain strong convergence of 
the density approximation. The following argument uses only the continuity 
approximation and is very similar to the corresponding argument in [14, 15, 
16]. We include it here for the sake of completeness. 

Lemma 10.1. Let (gh,Uh) be the numerical solution obtain through Defini- 
tion 3.1 and (3.7). The density approximation converges a.e 

Qh —> g a.e on (0, T) x Q. 

Proof. From Lemma 4.1 with B(g) = glog g, we have the inequality 

/ Qh log Qh dx (t) - / g° h log g° h dx < - g h div u h dxdt. 

Jn Jn Jo Jn 

From Lemma 2.3, we know that the limit (g, u) is a renormalized solution 
to the continuity equation. In particular, the following identity holds 

/ g\ogg dx (t) — / £olog£o dx = — / p div it dxdt. 
Jn Jn Jo Jn 

Note that there is no problems with integrability of g div u. Indeed, since g £ 
L°°(0,T;L^(Q.)) with 7 > 3 and divu G L 2 (0, T; L 2 (S1)), Holders inequality 
provides gdivuG L 2 (0, T; Ls (0)). 

By subtracting the two previous identities, we obtain 

/ gh log Qh - g log g dx (t) 
Jn 

< / g div u - g h div h u h dxdt + / g° h log g° h - g log £ dx 
Jo Jn Jn 

= / (f> 2 (g div u — gh div^ Uh) dxdt 
Jo Jn 

+ / / (1 - 4> 2 ){g div u - g h div h Uh) dxdt 
Jo Jn 

+ I Qh log Qh ~ Qo log go dx. 
Jn 

From Proposition 9.1, we have that 

lim / / (f) 2 (g div u — gh div Uh) dxdt 
h-^oJ J n 



h^0 



Jn 



lim / / <j) (gp(g) - ghP(gh)) dxdt < 0, 



where the last inequality follows from the convexity of p(g). Hence, this, 
together with the strong convergence of the initial density, allow us to con- 
clude 



/ g log g - g log g dx (t) < / / (1 
Jn Jo Jn 



4> 2 )Q dxdt, 
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where Q = gdivu — q&ivu € L 2 (0, T; L$ (f2)). Hence, we see that 



0< / glog g - glog dx (t) < C \\l - (f> 2 \ 



L 6 (n) ' 



where the first inequality follows by convexity of z H > z log z. Now, for any 
e > 0, we can choose <fi G W '°°(0) such that 



< q log £> — log <ix (i) < e. 



As a consequence, £>log £> — £?log q = a.e and hence Qh ^ Q a.e in (0, T) x 
O. □ 

10.1. Proof of the main result (Theorem 3.5). To conclude the main 
result, the only remaining part is to prove that (q, u) is a weak solution 
of the momentum equation (1.2) and that it satisfies the energy inequality 
(2.1). The other statements in Theorem 3.5 are all covered by (7.1), Lemma 
7.1, Lemma 7.2, and Lemma 10.1. Since we now know that Qh —> Q a.e, we 
have in particular 

p{q) = p(q) a - e i n (0)^) x 

By applying this information in (7.7), we immediately see that (q,u) is a 
weak-solution of the momentum equation. 

By passing to the limit h — > in the numerical energy inequality (4.1) 
(Lemma 4.2), using convexity, we discover that the limit (q,u) satisfies the 
energy inequality (2.1). 



Appendix A. Existence of a numerical solution 

Since the numerical method in Definition 3.1 is nonlinear and implicit it 
is not trivial that it is actually well-defined (i.e admits a solution). In addi- 
tion, the discretization of the momentum transport is posed using element 
averages of the velocity. As we will see the only part of the discretization 
that provides sufficient number of equations to determine all the degrees of 
freedom of Uh is the discretization of the diffusion operator. Hence, in it's 
present form, our discretization is not suitable for the Euler equations. 

The purpose of this section, is to prove the following the existence result 
which we have relied on in our analysis. 

Proposition 3.3. For each fixed h > 0, there exists a solution 

( e k h ,u k h )eQ h (n)xv h (n), e k h (-)>o, k = i,...,M, 

to the numerical method posed in Definition 3.1. 

To prove this result, we shall use a topological degree argument. The 
argument is strongly inspired by a very similar argument in the paper [12]. 
We will argue the existence of solutions to the following finite element map. 

Definition A.l. Let the finite element map 

H : Q+(0) x V h (n) x [0, 1] i — ^ Q h (n) x V h (n) be given by 



H(g h ,u h ,a) = (f h (a),g h (a)), 
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where (fh(a), gh(a)) are obtained through the mappings: 

-_A 

At 



Qh — Q k ~ l 
fh(a)q h dx = I — q h dx 



n 

— a 



^^ U P(H Wr dS ( x ) (A.l) 
+ o/ l 1 - £ ^| r [ a ] r W r dS(x), 



for all qh £ Qh(&) and 



g h (a)v h dx = I ^-^ — — ^— v h dx + V h u h V h v h dx 
n AJ Jo. 



a J \Jp(gu ® 2) [vj r dS(x) 



(A.2) 



for all G V/^fi). 



Observe that a solution of H(gh,Uh, 1) = (0,0) is a solution to our nu- 
merical method as posed in Definition 3.1. 

Before proceeding, let us make clear what we mean by topological degree 
in the present finite element context and denote by d(F, 0, y) the Z-valued 
(Brouwer) degree of a continuous function F : O — > M. M at a point y 6 
R M \F(dO) relative to an open and bounded set O C R M . 

Definition A.2. Let Sh be a finite element space, || • || be a norm on this 
space, and introduce the bounded set 

Sh = {Qh 6 S h ; \\q h \\ < C} , 

where C > is a constant. Let {o~i}f£ 1 be a basis such that spanjo-j}^ = Sh 
and define the operator Ilg : Sh — > M M by 

M 

H-B<ih = (qi,Q2,---,qm), qh = ^2qi°i- 

i=l 

The degree ds h (F, Sh,qh) of a continuous mapping F : Sh — > Sh at G 
Sh\F(dSh) relative to i?/, is defined as 

d Sh (F,S h ,q h ) = d(u B F(U B 1 ),U B Sh,n B q h 

The next lemma is a consequence of some basic properties of the degree, 
cf. [7]. 

Lemma A. 3. Fix a finite element space Sh, and let ds h (F, Sh, qh) be the 
associated degree of Definition A.2. The following properties hold: 

(1) ds h (F, Sh,qh) does not depend on the choice of basis for Sh- 

(2) d Sh (ld,S h ,q h ) = l. 
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(3) ds h (H(-,a),Sh,qh(,&)) is independent of a £ J := [0,1] for H:Sh X 
J Sh continuous, qh : J — >■ Sh continuous, and qh(ct) ^ H(dSh, a) 
Va€[0,l]. 

(4) d Sh (F,S h ,q h ) ^Q^F- l (q h ) ^ 0. 

To prove Proposition 3.3, we shall apply Lemma A. 3 with q^ = and 
mapping H given by Definition A.l. Let us first prove that our mapping H 
satisfies (3) in Lemma A. 3. 

Lemma A.4. Let H : Q^(fy x V h (tt) x [0, 1] i-> Q h (£l) x V h (Vt) be the finite 
element mapping of Definition A.l. There is a subset Sh C x V^(O) 

for which H : Sh X J — > Sh is continuous and the zero solution (0, 0) 
H(dS,a) for all a £ [0,1]. 

Proof. For any subset S C Q^(Q) x Vh(£l) bounded independently of a, 
the corresponding mapping H(Sh, a, 0) is clearly continuous. This follows 
directly from (A.l) and (A. 2) using the equivalence of finite dimensional 
norms. The more involved part is to determine a subset Sh for which (0, 0) 
H(dS, a) independently of a. 

Now, let us for the moment assert the existence of of (g, u) satisfying 

H(g h ,u h ,a) = (0,0). 

Then, from Lemma 3.4, we have that Qh > and moreover (A.l) yields 

/ Qh dx = I Q k h ~ l dx. 
Jn Jn 

Consequently, we can conclude that 

< C t , (A.3) 

independently of a. 

To derive a bound on the velocity Uh, we can repeat the steps in the proof 
of Proposition 4.2 (the energy estimate) while keeping track of a to obtain 

2 



I QjAu^ + a dx + At f ]Vuh] 2 dx 

Jn 2 7-1 Jn 

Jn 1 7—1 



< / 6h Z h 1 + -^-tpU- 1 ) dx < C, 



2 ' 7 -T^~ 

where C is independent of a. Together with (A.3), this allow us to conclude 

||£/ilk°°(fi) + ll n /i||L°°(n) < Cf 
We can now define the subspace 

Sh = {(Qh,Uh) e Ql x V h ; 1 1 g h 1 1 L oo (q) + H^llioofn) < c t } , 

which by definition has the property that (0, 0) H(dS, a) for all a £ [0, 1]. 
This concludes the proof. □ 
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Lemma A. 5. Let S be the subspace obtained by the previous lemma. Then, 
the topological degree of H(Sh,0) at qh = is non-zero: 

d Sh (H(-,0),S h ,0)^0. (A.4) 
As a consequence, there exists (gh,Uh) £ Sh such that 

H( 6h ,u h ,l) = (0,0), 
and hence Proposition 3.3 holds true. 

Proof. First, we note that proving (A.4) is equivalent to proving the exis- 
tence of (gh,u h ) € xV h satisfying, for all (qh,Vh) GQfcX V h , 



/ g h q h dx = g k h V dx, 
Jn Jn 

/ QhUhVh dx + At / V h u h V h v h dx = Q k h ~ 1 u k h ~ 1 v h dx. 
Jn Jn Jn 



(A-5) 



The first equation has the solution Qh = q\ 1 ■ Setting this into the second 
equation in (A. 5), we see that the resulting linear system is a sum of a posi- 
tive matrix Q k ~ 1 UhVf l and a symmetric positive definite matrix AtV hUhS 'hVh- 
Since the Laplace problem with the Crouzeix-Raviart element space and 
dirichlet conditions is well-defined, there is no problems with concluding the 
existence of Uh satisfying the second equation in (A. 5). 

□ 
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